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PREFACE 


The  author  has  been  involved  in  the  simulation  of  guided  projectiles  for  many 
years.  Different  investigators  adopt  different  coordinate  frames,  such  as  body- 
fixed,  plane-fixed  and  aeroballistic  (zero  P).  Some  use  the  Euler  angle 
representation  to  deal  with  rotations  and  some  use  quaternions.  Sources  which 
explain  the  significance  and  advantages  and  disadvantages  of  these  various 
approaches  are  not  readily  available.  It  is  difficult  to  find  derivations  and  there 
is  a  lack  of  advice  on  incorporating  these  methods  into  computer  simulations. 
The  author  became  especially  frustrated  when  he  attempted  to  collect  the 
equations  to  convert  an  existing  six  degree  of  freedom  (6  DOF)  simulation  from 
the  Euler  angle  to  the  quaternion  representation.  Several  sources  for  the  needed 
equations  were  found  but  no  two  agreed  exactly.  Since  little  in  the  way  of 
derivations  were  provided,  it  was  not  trivial  to  verify  the  equations  or  reconcile 
the  discrepancies.  This  document  resulted  from  the  author’s  attempt  to  make 
some  sense  of  this  confusion. 

The  first  chapter  contains  an  overview  of  the  problem.  In  the  second  chapter,  a 
brief  review  of  the  bare  minimum  of  matrix  algebra  is  provided  to  remind  the 
reader  of  some  of  the  important  properties  of  orthogonal  transformations.  The 
third  chapter  develops  the  Euler  angle  formalism  with  an  introduction  to  the 
difference  between  body-fixed,  plane-fixed  and  aeroballistic  coordinates.  The 
quaternion  algebra  is  developed  in  chapter  4.  This  is  an  extensive  subject.  Only 
enough  of  the  formalism  was  developed  to  provide  understanding  of  quaternions 
and  introduce  the  tools  needed  for  this  document.  In  chapter  5  the  rigid  body 
equations  of  motion  are  developed  for  the  three  coordinate  frames  discussed,  in 
both  the  Euler  angle  and  quaternion  representations.  The  discussion  of  the 
distinction  between  body-fixed,  plane-fixed  and  aeroballistic  coordinates  is 
distributed  throughout  chapters  3  to  5.  In  chapter  6,  the  integration  of  the 
equations  of  motion  is  discussed.  Discussions  of  the  treatment  of  Coriolis  and 
centripetal  corrections  in  a  flat  earth  model,  gravity  for  a  non-flat  earth,  and  time 
varying  mass  and  moment  of  inertia  have  been  included  in  this  report.  These 
topics  will  be  treated  in  a  future  report.  The  appendix  provides  a  summary  of 
the  algorithms  needed  for  implementing  these  results  in  a  6  DOF  simulation. 

The  author  wishes  to  thank  Dr.  Richard  Haddad  of  Polytechnic  University  of 
New  York,  Messrs  Romel  Campbell  and  John  Grau  of  ARDEC,  Dover,  NJ,  and 
Mr.  Thomas  Harkins  of  ARL,  Aberdeen,  MD,_  for  valuable  discussions  and 
.suggestions.  He  wishes  to  thank  Mr  Sung  Chung  for  checking  the  mathematical 
derivations. 
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1  INTRODUCTION 


When  developing  simulations  of  aircraft,  missiles  or  gun-launched  projectiles, 
investigators  require  a  coordinate  frame  in  which  to  follow  the  motion. 
Newton’s  laws  require  an  inertial  (unaccelerated)  frame.  The  earth  is  a 
convenient  reference  frame  but  is  not  inertial  since  the  earth  rotates.  The  earth 
may  nontheless  be  used,  with  Coriolis  and  centripetal  accelerations  included  to 
account  for  the  earth’s  rotation. 

However,  the  projectile  is  both  translating  and  rotating.  Thus  it  is  convenient  to 
express  the  equations  of  motion  of  the  projectile,  missile  or  aircraft  in 
coordinates  that  move  along  with  it  in  some  way.  The  obvious  choice  is  body- 
fixed  coordinates.  These  coordinates  are  attached  to  the  projectile  or  aircraft  and 
roll,  pitch  and  yaw  with  it.  The  reader  familiar  with  gimbals  or  gyroscopes  will 
recognize  that  these  Euler  angles  of  roll,  pitch  and  yaw  are  equivalent  to  gimbal 
angles.  In  the  case  of  a  guided  projectile,  the  seeker,  rate  sensor, 
accelerometers,  and  control  mechanisms  whether  aerodynamic  or  reaction  control 
all  operate  in  and  are  easiest  to  describe  in  body-fixed  coordinates. 

Sometimes  non  rolling  coordinates  are  desirable.  It  is  difficult  to  interpret  results 
of  a  simulation  when  the  point  of  view  is  rolling,  as  they  are  with  body-fixed 
coordinates.  In  addition,  spin-stabilized  gun  fired  projectiles  rotate  at  hundreds 
of  revolutions  per  second.  Computer  run  times  for  such  projectiles  using  body- 
fixed  coordinates  become  intolerably  long.  This  difficulty  arises  because  the 
integration  time  step  must  become  extremely  small  in  order  to  keep  the  angle  of 
roll  small  during  the  integration  time  step.  If  this  is  not  done,  gravity  is  smeared 
over  the  angular  motion  that  occurs  during  the  integration  time  step  because  of 
the  high  roll  rate,  giving  incorrect  results. 

Some  type  of  non-rolling  coordinate  system  is  used  to  deal  with  this  problem. 
One  solution  is  to  set  the  x  component  of  the  coordinate  frame  angular  velocity 
to  zero.  Another  is  to  set  the  Euler  roll  angle  to  zero.  These  two  approaches  are 
not  identical,  as  we  shall  see  in  subsequent  chapters.  We  shall  see  that  the 
difference  arises  from  the  fact  that  the  components  of  the  angular  velocity  form 
an  orthogonal  set  whereas  the  three  Euler  angles  do  not  have  a  mutually 
orthogonal  set  of  rotation  axes. 

Choosing  the  roll  Euler  angle  to  be  zero  eliminates  the  horizontal  component  of 
gravity  in  a  flat  earth  model  entirely  since  we  shall  see  that  the  y-axis  is 
constrained  to  move  in  a  horizontal  plane.  This  makes  the  numerical  integration 
insensitive  to  the  roll  rate.  However,  it  is  still  sensitive  to  the  pitch  and  yaw 
rates.  This  approach  is  typically  selected  when  modeling  an  unguided  stage  of  a 
spin  stabilized  projectile.  This  type  of  frame  is  called  plane-fixed. 

Choosing  the  x  component  of  the  coordinate  frame  angular  velocity  to  be  zero 
yields  aeroballistic  coordinates.  This  choice  does  not  completely  eliminate  the  y 
component  of  gravity  but  sensitivity  to  the  effects  of  roll  is  greatly  reduced.  Its 
chief  value  is  the  simplification  of  the  equations  of  motion.  Coupling  terms 
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involving  the  x  component  of  the  frame  angular  velocity  disappear  from  the 
equations  of  motion.  If  further  simplifications  are  made  based  on  symmetry  and 
linearity  of  the  aerodynamics,  it  is  possible  to  obtain  closed  form  solutions  to  the 
equations  of  motion1. 

With  any  of  these  frames,  it  is  necessary  to  regenerate  the  frame  as  the  projectile 
moves.  Thus  the  frame  itself  has  equations  of  motion.  We  shall  see  that  the 
rotation  matrix  that  transforms  the  vectors  from  the  moving  frame  to  the  inertial 
(earth)  frame  can  be  expressed  in  terms  of  either  three  Euler  angles  or  four 
quaternions.  The  equations  of  motion  for  the  Euler  angles  and  for  the 
quaternions  are  derived  so  that  they  may  be  integrated  to  obtain  the  new  frame 
and  update  the  projectile  equations  of  motion.  Only  two  angles  are  required  to 
describe  the  rotation  of  a  rigid  body  so  not  all  the  Euler  angles  or  the 
quaternions  are  linearly  independent.  Constraints  such  as  normalization 
conditions  therefore  exist  and  will  be  derived. 

The  advantage  of  Euler  angles  over  quaternions  is  their  intuitiveness.  Roll,  pitch 
and  yaw  are  a  natural  way  for  a  pilot  to  describe  or  visualize  the  angular  motion 
of  an  aircraft.  The  Euler  angles  are  the  natural  variable  for  describing  a  seeker 
or  spinning  gyroscope  gimbal.  However  the  Euler  angle  algebra  is  somewhat 
messy  and  unsymmetrical,  so  errors  are  not  always  evident.  Furthermore,  the 
sine  and  cosine  of  the  three  Euler  angles  must  be  repeated  computed,  providing  a 
computational  burden  that  does  not  exist  with  quaternions.  Thus,  although 
quaternions  are  not  intuitive  in  the  sense  that  Euler  angles  are,  their  simplicity 
and  symmetric  form  make  derivations  much  simpler,  are  less  prone  to  mask 
errors  and  are  computationally  more  efficient.  No  trigonometric  functions  or 
transcendental  functions  need  to  be  evaluated.  The  most  complicated  quaternion 
arithmetic  requires  the  square  of  a  quaternion  or  the  product  of  two  quaternions. 
For  this  reason  quaternion  algebra  is  desirable  in  digital  autopilots  for  guided 
projectiles  because  it  alleviates  the  computational  burden.  Furthermore,  Euler 
angles  are  susceptible  to  singularities  that  can  be  avoided  by  using  the  quaternion 
formalism. 

The  details  of  the  Euler  angle  and  quaternion  formalism  required  to  develop  the 
6  DOF  equations  of  motion  for  a  rigid  body  in  the  three  coordinate  frames 
discussed  above  will  be  developed  in  subsequent  chapters. 


1  Vaughn,  Harold  R.,  "A  detailed  Development  of  the  Tricyclic  Theory,"  Sandia  Laboratories,  SC-M-67-2933, 
Albuquerque,  NM,  1968. 


2  REVIEW  OF  MATRIX  ALGEBRA 
FOR  ORTHOGONAL  TRANSFORMATIONS 


This  chapter  contains  a  brief  review  of  the  matrix  algebra  required  in  this 
document1,2.  An  n  by  m  matrix  A  is  an  array  of  elements  a..,  i  =  1  to  n,  j  =  1  to 
m,  of  i  rows  and  j  columns,  which  obeys  the  following  ’/aws  of  addition  and 
multiplication. 


to 

+ 

◄ 

II 

u 

c. .  =  a. .  +  b . . 

*J  ij  *J 

(2.1) 

C  =  A  B 

- 

c. .  =  y  a,,  b ,  . 

ij  ik  kj 

(2.2) 

j 


For  these  operations  to  be  meaningful,  certain  matching  restrictions  exist  on  the 
number  of  rows  and  columns.  For  addition.  A,  B  and  C  must  all  have  the  same 
number  of  rows  and  the  same  number  of  columns.  For  multiplication,  the 
number  of  columns  of  A  must  match  the  number  of  rows  of  B.  The  product  C 
has  the  same  number  of  rows  as  A  and  the  same  number  of  columns  as  B.  Such 
matrices  are  said  to  be  conformable  .If  the  number  of  rows  and  columns  are 
equal,  the  matrix  is  square.  A  vector  can  be  represented  by  an  n  by  1  column 
matrix. 

The  usual  algebraic  laws  hold  except  that  multiplication  is  not  generally 
commutative  and  the  multiplicative  inverse  does  not  always  exist  (see  below). 
When  we  refer  to  the  inverse  of  a  square  matrix,  we  generally  mean  the 
multiplicative  inverse.  The  inverse  of  a  square  matrix  A  is  denoted  by  A  -1  and  is 
defined  by 

A  A_1  =  A_1  A  =  1  (2.3) 

where  1  is  the  unit  matrix  (i.e.,  1  along  the  diagonal  and  zero  elsewhere).  This 
can  also  be  written  in  terms  of  the  elements  of  the  matrix 


y  a. 

^  u 


- 1 


y  a. . 

^  ‘J 


jk 


=  8. 


ik 


(2.4) 


j 


} 


1  Wylie,  C.  R.,  Jr,  "Advanced  Engineering  Mathematics,"  McGraw-Hill,  New  York,  1956. 

2 

Margenau,  Henry  and  George  Moseley  Murphy,  "The  Mathematics  of  Physics  and  Chemistry,"  Van  Nostrand,  New 
York,  1956. 


where  i  i  and  k  -  1,  2,  3  and  is  the  Kronecker  delta  function  which  is  unity 
when  the  subscripts  are  equal  and  zero  otherwise.  In  general,  the  inverse  of  a 
matrix  does  not  always  exist.  Generally,  a  necessary  and  sufficient  condition  for 
^matrix  to  have  an  inverse  is  that  it  must  be  square  and  the  determinant  of  the 
i  7*fA  a  matrix  is  called  non  singular.  Even  when  a 

“«r!x  Tnon  singular  fining  .he  inverse  is  generally  tedious.  This  will  no. 

“oncern  us  here  since  we  will  be  dealing  only  with  orthogonal  matrices  wtahate 
concern  us  nerc  m.  .  ,  We  shall  also  see  that  the  inverse  can  be  found 

quite*  trivially  From  this  point  onwards,  we  adopt  the  summation  convention 
Jhfch ^states  that  a  sum  is  implied  over  any  index  that  is  repeated  twice.  Thus 


a.,  b  =  ^ ik 

ij  jk  "  ij  J* 


(2.5) 


preserve  length  since  rotation  does  no^tretch  gn™jr^^°j^s^\phat 

^operelrt/derCm  r^sVeSion,  Confider  the  matrix  A  operating 
on  the  column  matrix  or  vector  v. 

,  a  (2.6) 

v'  =  A  v 


This  may  also  be  written  in  terms  of  the  elements 


v '  =  a.,  v. 
»j  J 


If  the  length  of  v  is  preserved  by  this  transformation  A,  then 

v' .  v'  =  v.v 


(2.7) 


(2.8) 


where  the  dot  product  is  defined  by 


2  2,2 

V.V  =  V,  +  V_  +  V 


Thus,  if  A  is  an  orthogonal,  length-preserving  transformation 


v  V;  =  v'.  v.  -  a.. v.  a.kvk  a..a.kv.vk 


it  i  i  >J  J 


This  is  possible  only  if 


a.  a  v  =  6 
ij  ik  j* 


But  this  is  the  definition  of  the  inverse  of  A 


aji  1  V“  hjk 


(2.9) 


(2.10) 


(2.11) 


(2.12) 


Thus,  for  an  orthogonal  matrix,  aJ1  «  V  The  inverse  of  an  orthogonal  matrix 
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is  obtained  by  interchanging  the  subscripts  or  the  rows  and  columns  of  the 
matrix.  This  is  equivalent  to  reflecting  the  matrix  about  the  diagonal.  This 
operation  is  called  the  transpose  and  will  be  denoted  by  a  superscript  T.  Thus 

A-=  AT  (213) 

defines  an  orthogonal  matrix 

We  define  the  trace  or  spur  of  a  matrix.  The  trace  is  simply  the  sum  of  all  the 
elements  along  the  diagonal.  Thus 

Tr  [  A  ]  =  Tr  [  AT  ]  =  2an  (2‘14) 


The  trace  is  obviously  invariant  to  the  transpose  operation  since  the  diagonal 
elements  are  unchanged  under  transposition. 


Finally,  we  show  that  the  transpose  of  a  product  of  two  matrices  is  the  product  of 
the  transposes  of  the  individaul  matrices,  but  in  reverse  order. 

[  A  B  ]T  =  BT  AT  (2-15) 


The  proof  is  as  follows: 


[  A  B  ] 


T 

ik 


B  . 
q» 


XT  T  T 

B  A  =  B  A  =  [  B  A  ] 
qi  kq  iq  qk 


ik 


(2.16) 


1  It  can  be  shown  that  the  determinant  is  invariant  to  the  transpose  operation.  Therefore  from  (2.3)  and  (2.13)  we  can 
conclude  that  the  square  of  the  determinant  of  an  orthogonal  matrix  is  1.  Thus  the  determinant  must  be  ±1.  The 
negative  determinant  is  associated  with  reflections,  which  obviously  also  preserves  length.  The  positive  deter¬ 
minant  is  associated  with  rotations. 
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3  EULER  ANGLES 


3.1  ROTATING  COORDINATE  FRAMES 

We  use  a  right  hand  coordinate  system  with  x  positive  forward,  y  positive  to  the  right  and 
z  positive  down,  as  indicated  in  Figure  1 .  Each  Euler  angle  also  obeys  a  right  hand  rule 
with  respect  to  its  axis.  The  roll  <p  occurs  about  the  x  axis,  the  pitch  0  occurs  about  the 
y  axis,  and  the  yaw  y  occurs  about  the  z  axis.  Thus  the  roll  is  positive  clockwise  looking 
forward  from  the  rear,  pitch  is  positive  upward  (even  though  z  is  positive  down),  and  the 
yaw  is  positive  looking  forward.  By  convention,  the  rotational  transformation  from  iner¬ 
tial  to  rotating  axes  consists  of  a  yaw  through  angle  \|/,  followed  by  a  pitch  through  angle 
0  and  finally  a  roll  through  angle  <p.  The  order  of  these  rotations  is  important  since  ro¬ 
tations  do  not  commute.  The  component  rotations  are  shown  in  Figure  1.  The  original 
inertial  axes  are  indicated  by  x,  y,  and  z.  The  primes  indicate  intermediate  axes  of 
subsequent  rotations.  The  final  rotating  axes  are  indicated  by  x”,  y”\  and  z’”. 

These  component  rotations  may  be  expressed  as 


cos(t|0 

sin(t|i) 

0 

=  —  sin(tJO 

cos(t)t) 

0 

(3.1) 

0 

0 

1 

cos(O) 

0 

—  sin(0) 

R 

-  0 

1 

0 

(3.2) 

e 

sin(0) 

0 

cos(0) 

1 

0 

0 

** 

=  0 

cos(<|>) 

sin(4») 

(3.3) 

0 

—  sin(4>) 

cos(4>) 

1  Hakelock,  John  H-,  "Automatic  Control  of  Aircraft  and  Missiles,"  John  Wiley  and  Sons,  New  York,  1965. 

2  Ptvin  Bernard,  "Dynamics  of  Flight  -  Stability  and  Control,"  John  Wiley  and  Sons,  New  York,  1965. 


3  Etlrin,  Bernard,  "Dynamics  of  Atmospheric  Flight,"  John  Wiley  and  Sons,  New  York,  1972. 
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Figure  2.  Component  Euler  Angle  Rotations 


0  :£  4>  <  2ir  , 


—  IT 


(3.4) 


ir 

<  8  <  ,  —  TT  <  l|/  <  TT 

2  2 

It  is  conventional  to  think  of  a  rotation  from  the  non-moving  (e.g.,  inertial) 
coordinates  to  the  moving  body-fixed  coordinates  in  terms  of  a  yaw  followed  by 
a  pitch  and  finally  a  roll.  However,  the  more  commonly  used  transformation  in 
simulations  will  be  the  inverse  rotation  from  body-fixed  to  non-moving 
coordinates.  To  find  this  transformation,  we  must  reverse  both  the  order  of  the 
component  subrotations  and  reverse  the  algebraic  sign  of  the  roll,  pitch  and  yaw 
angles.  Thus  the  matrix  rotation  that  transforms  a  vector  from  the  rotating  axes 
to  the  inertial  axes  is  obtained  by  evaluating  the  matrix  product 


cos0  cosi|* 


—  cost))  sin»|>  sin4>  sintj* 

4-  sin<J>  sin0  cos»|>  +cos4>  sin8  cosi|* 


cos0  sim|» 


cost})  cost}*  —  sin<}>  cost}; 

+  sin4>  sin0  sirnj*  +cos<}>  sin0  sinij* 


(3.5) 


—  sin0 


sin<}>  cos0 


cos<}>  COS0 


Since  this  matrix  is  a  rotation  and  therefore  orthogonal,  the  inverse 
transformation  from  inertial  to  rotating  coordinates  is  obtained  by  taking  the 
transpose,  i.e.,  interchanging  rows  and  columns.  Later  on  we  shall  distinguish 
between  various  types  of  moving  coordinate  systems:  body-fixed,  plane-fixed, 
and  aeroballistic  coordinates. 

Another  interesting  set  of  expressions  that  is  needed  is  the  components  of  the 
angular  velocity  ft  of  the  coordinate  axes  in  terms  of  the  Euler  angles  and  their 
derivatives.  This  will  be  useful  in  the  development  of  the  equations  of  motion. 
An  uncritical  guess  might  be  that  the  components  of  ft  are  nothing  more  than  the 
derivatives  of  .the  three  Euler  angles.  While  it  is  true  that  ft  is  the  vector  sum  of 
the  vectors  associated  with  these  Euler  angle  rates,  these  rate  vectors  are  not 
mutually  orthogonal  and  so  can  not  be  the  components.  However,  they  can  be 
expressed  in  such  a  way  that  they  can  be  resolved  into  orthogonal  components  in 
the  moving  frame.  In  general,  the  rotation  associated  with  ft  can  be  considered 
as  consisting  of  three  successive  non-orthogonal  rotations  with  angular  velocities. 
See  Figure  2.  In  the  following,  a  bar  denotes  a  vector  and  a  hat  denotes  a  unit 
vector. 
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(3.6) 


XI  =  XI  +  XI  +  XI 

4/  8  <}> 

=  4*  f  +  e  >•'  +  <t>  x" 

We  want  to  resolve  this  in  moving-fixed  coordinates,  where  we  can  write 

a  =  a  „  x"  +  a  y'"  +  a  £"'  (3.7) 

x  y  z  7 

The  vector  quantities  on  the  right  of  (3.6)  are  not  an  orthogonal  set  whereas 
those  in  (3.7)  are  orthogonal.  The  unit  vectors  x" ,  y and  £"'  are 
orthogonal.  These  may  be  resolved  as  follows.  The  vector  a^  resolved  in  the 
orthogonal  coordinates  in  the  moving  frame  is  obtained  by  applying  the  partial 
Euler  rotation: 


(o  ^ 

1. 0 

' 

1 

(n  ^ 

=  !  0 

6  e  1 

y"‘ 

i 

(n  ^ 

\r*)  . 

z" ' 

i  * 

I  -i|i  .«n(9) 

=  I  »j/  sin  (<f>)  cos  (0) 
I  co.y(<t>)  co.v(0) 

Likewise, 


K) 


(3.8) 


(3.9) 
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0 


+  0  cos( 4* ) 
—  6  sin(4>) 


No  transformation  is  needed  for  the  components  of  since  it  is  parallel  to  the 
x"-axis.  Thus 


K)  -  * 

X  ' 

K) ...  -  0 


(3.10) 


(3.11) 


(3.12) 


We  can  add  vectorially  (3.8)  to  (3.10)  to  obtain  an  explicit  representation  of 
(3.6)  resolved  in  the  moving  frame.  Comparing  this  to  (3.7)  gives  the  result  we 
seek,  viz. 


—  —  5t/i(0)  +  <f> 

(3.13) 

ft 

y 

=  t}>  sm(< |>)  cos  (6)  +  8  co^(<})) 

(3.14) 

ft 

2 

=  co5,(4>)  cos( 9)  —  0  sin(<j>) 

(3.15) 

These  can  be  inverted  for  the  derivatives  of  the  Euler  angles  by  the  usual 
algebraic  techniques  to  get  (dropping  the  primes) 

[fty  sin  (<j>)  +  ft?  cos(c|>)  j 

i  =  (3.16) 

cos (9) 


it 


(3.17) 


4>  =  +  [n^  .?in(4»)  +  flz  co.v(<J>)j  tan(d) 


0  =  [ly  cos(i J>)  -  ft2  sin(<$>)  (3.18) 

In  a  computer  simulation,  equations  (3.16)  through  (3.18)  are  numerically 
integrated  to  update  the  three  Euler  angles.  These  are  used  in  turn  to  update  the 
rotation  matrix  (3.5).  Note  that  the  above  expressions  have  a  singularity  at  0  = 
-ir/2.  Euler  angles  can  be  chosen  so  that  the  singularity  occurs  elsewhere,  but 
there  is  always  a  singularity.  In  simulations  using  Euler  angles,  care  must  be 
taken  that  the  vicinity  of  the  singularity  is  avoided. 


By  inspection,  (3.16)  through  (3.18)  can  be  put  into  matrix  form. 


(3.19) 


(3.20) 


Expressions  (3.19)  and  (3.20)  are  not  symmetrical  or  elegant.  Note  that  the 
matrices  in  (3.19)  and  (3.20)  must  be  inverses  of  one  another.  This  can  be 
verified  by  taking  the  product  of  these  two  matrices  and  verifying  that  the  unit 
identity  matrix  results. 


1  0  -sin(0)  <t> 

ftv  =  0  cos(d>)  sin(d>)cos(0)  0 

SI.  0  -sin(<)>)  cos(d>)cos(0)  vp 

Likewise,  (3.16)  through  (3.18)  may  be  written  from  inspection  as 


4> 

1 

sin(6)tan(9) 

cos(<j>)tan(0) 

SI 

X 

0 

= 

0 

cos(<f>) 

-sin(<j>) 

4- 

0 

sin(6)/cos(9) 

cos(<t>)/cos(0) 

a 
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The  choice  for  the  coordinate  frame  angular  velocity  components  ft  ,  ft  ,  and  ftz 
will  depend  on  the  application.  Body-fixed  coordinates  are  appropriate  for 
simulating  a  guided  projectile,  rocket  or  aircraft.  In  a  guided  projectile  the 
seeker,  sensors  and  control  mechanisms  are  fixed  to  the  body.  For  this  reason  it 
is  easiest  and  simplest  to  describe  these  subsystems  in  a  coordinate  frame  fixed  to 
the  projectile  body.  The  coordinate  frame  angular  velocity  components  are  then 
equal  to  the  analogous  components  of  the  body  angular  velocity,  which  are 

conventionally  denoted  by  P,  Q  and  R.  Thus  in  equations  (3.13)  through  (3.20), 
we  make  the  substitutions 

ft  =  P 

X 

ft  =  Q  (3.21) 

y 

ft  =  R 

Z 

This  is  the  usual  choice  made  for  6  DOF  simulations  of  guided  projectiles  and 
missiles. 


For  unguided  projectiles,  a  non-rolling  coordinate  frame  is  often  preferred.  Such 
a  frame  pitches  and  yaws  with  the  projectile  but  does  not  roll  with  it.  A  non¬ 
rolling  frame  might  be  defined  by  letting  the  x  component  of  the  frame  angular 
velocity,  ft  ,  vanish  or  by  letting  the  time  derivative  of  the  roll  Euler  angle,  <f>, 
vanish.  From  (3.13)  we  see  that  these  two  approaches  are  not  identical.  The 
coordinates  obtained  by  the  first  approach  are  called  aeroballistic  coordinates  and 
the  latter  plane-fixed  coordinates.  The  advantage  of  the  former  is  he 
simplification  of  the  equations  of  motion  since  the  coupling  terms  involving  ft^ 
drop  out,  as  we  shall  see  in  Chapter  5.  This  approach  is  often  taken  with  analytic 
or  closed-form  solutions  of  the  equations  of  motion.  Equation  (3.21)  would  be 
modified  by  letting  ft^  =  0. 

The  plane-fixed  approach  is  often  used  for  6  DOF  computer  simulations  of  spin 
stabilized  projectiles.  Spin  stabilized  projectiles  have  typical  spin  rates  of 
hundreds  of  revolutions  per  second.  Using  a  body-fixed  representation  in  a 
computer  simulation  of  such  a  projectile  requires  an  extremely  small  integration 
time  step  and,  consequently,  inordinately  long  computer  run  times.  The  time  step 
must  be  small  so  that  the  projectile  roll  is  not  appreciable  during  the  time  step. 
Otherwise  the  effect  of  gravity  is  smeared  across  this  angle.  While  aeroballistic 
coordinates  will  help,  a  more  useful  solution  is  to  require  d$/dt  =  cf>  =  0  for  the 
coordinate  frame.  We  shall  see  in  Chapter  5  that  the  y  component  of  gravity  is 
rigorously  eliminated  in  the  fixed  plane  approach.  This  eliminates  sensitivity  of 
the  integration  to  the  projectile  roll  rate,  though  a  similar  sensitivity  is  still 
present  for  the  much  slower  pitch  and  yaw  rates.  This  approach  can  speed  up 
simulation  run  time  by  orders  of  magnitude. 

All  three  approaches  will  be  discussed  further  when  developing  the  equations  of 
motion  in  Chapter  5.  The  plane-fixed  coordinates  are  derived  in  the  following 
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section  from  the  physical  view  point  of  constraining  one  of  the  axes  to  move  in  a 
single  plane. 


3.2  PLANE-FIXED  COORDINATES 


Plane-fixed  coordinates  pitch  and  yaw  with  the  body  but  do  not  roll  with  the 
body.  Hence  we  define  <{>  =  0  and  d$/dt  =  0.  More  precisely,  one  axis  is 
constrained  to  always  remain  in  one  plane,  though  it  can  rotate  in  that  plane.  For 
example,  the  z-axis  could  be  constrained  to  the  vertical  plane  (original  x-z 
plane).  This  can  be  achieved  in  an  inertial  to  moving  (i.e.,  plane-fixed)  frame 
transformation  consisting  in  a  pitch  about  the  original  y-axis  (which  keeps  the  z- 
axis  in  the  original  pitch  plane,  which  is  vertical)  followed  by  a  yaw  about  the 
new  z-axis,  which  leaves  the  z-axis  unchanged  and  therefore  still  in  the  vertical 
plane. 

Alternatively,  the  y-axis  can  be  constrained  to  the  horizontal  plane  (original  x-y 
plane).  This  is  done  by  a  yaw  about  the  z-axis  followed  by  a  pitch  about  the  y- 
axis.  We  can  construct  the  rotation  matrix  as  before,  using  (3.1)  and  (3.2). 
However,  unlike  before,  the  above  recipe  involves  inverse  transformations  from 
the  inertial  to  the  moving  frame  rather  than  from  moving  frame  to  inertial,  as 
was  the  case  in  the  derivation  of  (3.5)  for  the  body-fixed  frame.  Thus  the  above 
transformation  is  comprised  of  the  inverses  of  the  matrix  operators  previously 
used.  Thus,  recalling  (2.15)  and  that  the  inverse  of  an  orthogonal  matrix  is  its 
transpose. 


T  1  =  R  1  R  1 

0  il* 


T  T  T  T 

R  R  =  [  R  R  1  =  T 

e  >|)  1  >!>  e  J 


or 


T  =  R  ,  R  = 

i>/  e 


cos0  cosi|<  —  sini{i 

cos0  sintjj  cosv}» 

—  sin0  0 


sin0  costy 
sin0  sim|i 
COS0 


(3.22) 
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This  is  the  plane-fixed  analog  of  (3.5).  It  is  not  surprising  that  this  is  equivalent 
to  making  <£  vanish  in  (3.5).  Likewise  (3.6)  becomes 

n  =  a  +  n  =  z  +  e  y'  (3.23) 

0 

Equation  (3.8)  becomes 


° 

I 

0  I 
4»  ! 

(3.24) 

~i|/  sin(Q)  | 

6  I 

+  »(/  cos  (0)  I 


ft,  = 

4 


K) 

K) 


=  R. 


Likewise,  (3.9)  becomes 


1 

1 

(n  ^  1 

re)  I 

1 

0 

1 

X* '  j 

1 

| 

(n  )  1 

re)  1 

e 

I 

y'"  \ 

1 

1 

o  ! 

1 

0 

Equations  (3. 13)  through  (3.15)  become 
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ft  = 

z 

»Ji  cos(Q) 

(3.26) 

V  - 

e 

(3.27) 

n..  = 

—  ij;  sin(Q)  =  —  tan(0) 

(3.28) 

We  have  substituted  (3.26)  into  (3.28).  The  primes  on  the  subscripts  of  ft  in 
(3.26)  through  (3.28),  as  well  as  in  (3.13)  through  (3.15)  can  be  dropped  since 
the  axes  referred  to  are  orthogonal.  Inverting,  with  some  algebra  yields  the 
analogs  of  (3.16)  through  (3.18),  viz.. 

ft  —ft 

Z  X 

cos ( 9)  sin(8) 

(3.29) 

4> 

0 

(3.30) 

e 

ft 

V 

(3.31) 

The  analogs  of  (3.19)  and  (3.20)  are 

(lx  1  0  — sin(0)  0 

(3.32) 

n.  =  0  1  0  6 

ft,  0  0  cos(0)  ii» 

and 

0  1  0  tan(0) 

6=010  nv 

(3.33) 

4*  oo  i/cos(6)  n. 
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In  summary,  the  representation  of  the  rotation  matrix  T  is  unique  no  matter  how 
it  is  derived,  although  some  methods  may  not  always  work  because  of 
singularities.  The  rotation  matrix  T  can  be  viewed  as  an  operator  which  rotates  a 
vector  in  a  fixed  coordinate  system  or,  conversely,  as  a  rotation  of  the  coordinate 
system  while  the  vector  remains  fixed.  In  the  former  point  of  view,  the  vector 
has  the  same  length  but  its  components  are  changed  because  its  direction  in  space 
has  changed  due  to  the  vector’s  rotation.  In  the  latter  point  of  view,  the  vector 
has  the  same  length  and  direction  in  space  but  its  components  are  different 
because  of  the  rotation  of  the  coordinate  frame. 

If  we  take  the  latter  point  of  view,  we  can  see  that  the  rotation  matrix  is  just  the 
matrix  of  the  direction  cosines.  Let  i'  denote  the  three  mutually  orthogonal  unit 
basis  vectors  of  the  primed  (rotated)  coordinate  system  and  zAq  denote  the  three 
basis  vectors  of  the  unprimed  coordinate  frame.  T  will  be  the  elements  of  the 
transformation  from  the  unprimed  to  the  primed  frame.  Then 


=  T  i 
pq  q 


(3.34) 


Taking  the  dot  product  and  making  use  of  the  mutual  orthogonality  of  the  basis 
vectors,  we  can  write 


=  T  5  , 

pq  ql 


(3.35) 


Thus  the  elements  of  the  rotation  matrix  could  be  obtained  by  taking  the  dot 
products  of  the  unit  basis  vectors  of  the  unprimed  coordinate  frame  with  the 
basis  vectors  of  the  primed  frame. 
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4  QUATERNIONS 


We  will  develop  just  enough  of  the  algebra  of  quaternions  as  is  needed  for 
understanding  and  writing  six  degree-of-freedom  (6  DOF)  simulations. 
Quaternions  are  a  quadruplet  of  numbers  (strictly  speaking  operators)  that  can 
be  considered  to  be  a  generalization  of  complex  numbers.  Recall  that  the 
quantity 

i  =  V“i  (4.1) 

may  be  best  thought  of  as  a  rotation  operator.  Thus  i  is  a  90  degree 
counterclockwise  rotation  from  the  “real”  to  the  “imaginary”  axis.  See  Figure  3. 
The  "square"  of  i  is  two  successive  90  degree  rotations.  This  is  equivalent  to  a 
180  degree  rotation.  This  takes  us  to  the  negative  real  axis.  The  cube  of  i  is  a 
270  degree  counterclockwise  rotation  from  the  positive  real  to  the  negative 
imaginary  axis.  It  is  in  this  sense  that  i2  =  -1  and  i3  =  -i-.  The  4th  power  of  i  is 
just  a  360  degree  rotation  which  gets  us  back  to  the  real  axis. 

For  quaternions  we  define  three  such  quantities,  corresponding  to  rotations  about 
the  x,  y,  and  z  axes  respectively.  As  in  the  case  of  the  conventional  imaginary 
number  i,  the  operators  i,  j,  and  k  -can  be  interpreted  as  90  degree  rotations 
about  the  x,  y,  and  z  axes;  and  the  squares  correspond  to  a  180  degree  rotation 
about  the  appropriate  direction,  and  so  forth.  These  operators  are  sometimes 
call  hyperimaginary  numbers.  Just  as  the  conventional  complex  numbers  can  be 
used  to  provide'machinery  for  treating  rotations  in  a  plane,  we  might  expect  that 
three  “imaginary”  operators  i,  j,  and  k  might  be  used  to  treat  rotations  about 
three  axes,  i.e.,  in  three  dimensional  space.  Another  useful  property  of 
quaternions  is  that  it  permits  us  to  multiply  and  divide  vectors.  This  will  be  seen 
to  provide  a  much  simpler  mathematical  treatment  than  matrices. 

Recall  that  rotations  are  not  numbers  but  operators  and  do  not  commute.  Thus  ij 
does  not  equal  ji,  and  so  forth.  A  little  bit  of  thought  and  some  experimentation 
with  rotations  will  convince  the  reader  that  the  following  elementary 
relationships  hold. 


i  =  j 

=  k  =  -1 

(4.2) 

i  J  = 

k  =  -ji 

(4.3) 

j  k  = 

i  =  -k  j 

(4.4) 

k  i  — 

j  =  -/  k 

(4.5) 
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Figure  3.  Imaginary  Number  i  Interpreted  as  a  Rotation  Operator 
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We  take  a  quaternion  to  be  a  quadruple  Q  =  [  Qq  +  iQ  ^  +  jQ 2  +  kQ^].  Some 
tedious  algebra  will  verify  that  multiplying  out  two  quaternions,  taking  all 
possible  products  and  simplifying  using  (4.2)  through  (4.5),  gives  the  following 
result. 


W  Q  =  [  WQ+iW1  +  jW2+kW3  )  [  QQ+iQ1  +  jQ2+kQ3  ]  (4.6) 

=  [WoQo-W^.-WzQ.-W.Q.]  +  W  0[iQ  x  +  jQ  2  + kQ  3]  +  Q  0[iW  x  + jW  2+ kW  3) 

+  /[  W2Q.  -  W3Q2  ]  +  ;[  W3Q ,  -  WxQ3  ]  +  *[  W,Q2  -  W2Q1  ] 


Upon  inspection  of  (4.6),  we  see  that  the  first  line  after  the  last  equality  contains 
expressions  that  resemble  a  dot  product  and  a  cross  product.  This  suggests  an 
alternative  formulation.  We  can  treat  i,  j  and  k  not  as  rotation  operators  or 
hyperimaginary  numbers  but  as  an  orthogonal  set  of  unit  vectors  and  consider  a 
quaternion  formally  to  consist  of  a  scalar  part  and  a  vector  part.  This  lets  us 
write  the  quaternion  product  expressed  in  (4.6)  more  compactly. 

The  quaternion  product  is  equivalently  defined  as 

W  Q  =  [  WqQq  -  W.Q  ,WqQ  +  QqW  +  WxQ  ] '  (4.7) 

where 

W  «  [  WQ  ,  W  ]  and  Q  -  [  Q0  ,  Q  )  (4.8) 

Furthermore,  we  can  obtain  still  another  alternate  form  since  (4.6)  can  also  be 
rearranged  into 


WQ=  [  WoG0  “  Wi<2i  ~  W2Q2  "  W3Q3  ] 

+  i[WiQ0+W0Q1-W3Q2+W2Q3)  (4.9) 

+  /[  w2q0  +  w3q1  +  w0q2  -  w3q3  ] 

+  *[  -  WM,  +  WM -  +  WnQ,  ] 

L  3VJ0  2^1  1^2  0^3  J 

This  may  be  organized  into  the  matrix  form  as 
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W  Q  =  A  = 


*0 

+  w0 

-w2 

-W. 

Q0 

+  w} 

+  w0 

-W 3 

+  W2 

Qy 

A. 

+  W2 

+w3 

-Wy 

Q2 

*3 

+w3 

~w2 

+  w0 

e3 

+  Q0 

-Q: 

-q2 

“Cj 

^0 

+  Qy 

+  Go 

+  g3 

-e2 

+  G; 

-e3 

+  G0 

+  Gi 

W2 

+e3 

+g2 

~Qy 

+  G0 

W3 

(4.10) 


Compare  (4.7)  and  (4.10)  to  (4.6).  They  are  equivalent.  (4.7)  and  (4.10)  are 
more  compact  than  (4.6)  but  are  arbitrary  if  used  as  a  definition  for  quaternion 
multiplication,  as  some  authors  do.  On  the  other  hand,  (4.6)  flows  logically  and 
automatically  from  the  properties  (4.1)  to  (4.5),  and  gives  insight  into  the 
connection  with  rotations.  We  will  use  the  approach  that  is  the  most  convenient 
in  each  case. 


Some  useful  expressions  follow.  Define  the  conjugate 

Q*  -  t  Q0  , -Q  )  <4-n) 

Define  the  norm  or  absolute  value  squared 

IQ  I2  -  Q  Q*  =  [  q\  +  Q-Q  .  0  ]  =  Q*  Q  (4.12) 

Let  us  define  an  invese  and  verify  it  works. 

Q'1  ■  i  Q0  ,  ~Q  ]/(  q]  +  Q-Q  ) 

=  Q*  /  |  Q  |  (4.13). 

It  follows  that 
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(4.14) 


Q  Q'1  =  [  1  ,  o  ] 


If  the  norm  vanishes,  the  quaternion  is  said  to  be  singular  and  the  inverse  does 
not  exist.  It  is  easy  to  show  that  the  norm  of  a  product  equals  the  product  of  the 
norms.  The  inverse  of  a  product  is  the  product  of  the  inverses  in  reverse  order. 
The  conjugate  of  a  product  is  the  product  of  the  conjugates  in  reverse  order. 
Thus 


I Q2  Q2  1  “  IQjllQjl 

(4.15) 

[  Qj  Qj  1  1  -  Qj  '1  Qj'1 

(4.16) 

*  *  * 

[Q,  Qj]  -  Qj  Qj 

(4.17) 

In  general,  quaternion  arithmetic  will  be  familiar  except  for  non-commutativity 
of  multiplication.  Commutation  breaks  down  for  multiplication  because  of  the 
cross  product  term  in  (4.7).  Otherwise  all  the  other  usual  laws  are  obeyed. 
Quaternion  arithmetic  is  distributive  and  associative,  but  commutative  only  for 
addition.  Identity  elements  exist  for  both  addition  and  multiplication,  viz.  (0,0) 
and  (1,0).  Inverses  also  exist  for  addition  for  all  quaternions.  A  multiplicative 
inverse  exists  for  any  non-zero  quaternion.  See  (4.13).  The  rules  for 
differentiation  are  the  familiar  ones,  except  care  must  be  taken  because 

quaternions  do  not  commute.  As  an  example,  consider  the  derivative  of  the 
product  of  two  quaternions. 

[  Qj  Q2  y  =  Qj'  Q2  +  Qj  Q2'  *  Q2  +  Q2'  Qj  (4.18) 

Since  the  norm  of  a  product  equals  the  product  of  the  norms,  it  would  seem 

plausible  that  a  quaternion  of  unit  norm  could  be  useful  to  treat  rotations  since 

rotations  must  preserve  the  length  of  vectors.  We  shall  see  later  that  this 

conjecture  is  essentially  correct.  In  anticipation  of  a  quaternion  formalism  for 
rotation,  some  relations  for  unit  quaternions  will  now  be  provided. 


Consider  the  unit  quaternion  (sometimes  called  a  versor  or  Euler  quaternion) 

e  =  [  eQ  ,  e  ]  =  [  eQ  ,  ey  e2e3]  (4.19) 

where  |  e  |  2  =[1,0]  See  (4.12).  Thus 
2  2  2  2 

eQ  4-  e  +  e2  +  e  3  =  1  (4.20) 
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From  (4.15),  the  “length”  or  norm  of  a  quaternion  is  preserved  when  multiplied 
by  a  unit  quaternion. 

|eQ|-|Q|  (4-21) 


Since  the  norm  of  a  unit  quaternion  is  unity,  (4.13)  tells  us  that  the  inverse  of  a 
unit  quaternion  is  equal  to  its  conjugate,  i.e. 


* 


e 


(4.22) 


Also 


e  = 


. eo  ’  e\  e l  e3  . 


By  differentiating  e  e  =  1,  we  obtain 

.  -l  .-i 

e  e  =  —  e  e 


Applying  (4.16)  to  (4.24)  gives 

r .  -i  f1 

e  e  =  -  l  e  e  j 


.  - 1 

=  e  e 


- 1 


(4.23) 


(4.24) 


(4.25) 


Since  the  inverse  of  a  unit  quaternion  is  equal  to  its  conjugate,  the  inverses  in 
(4.25)  can  be  replaced  by  conjugation. 


r .  *  i 

e  e  =  -  [  e  e  J 

r  .*  ]' 

-  i“  \ 


(4.26) 


Vectors  can  be  treated  as  quaternions  with  a  zero  scalar  part.  Note  that 
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Some  authors  refer  to  a  quaternion  that  has  a  zero  scalar  part  as  a  pure 
quaternion.  Generally,  the  quaternion  product  of  two  pure  quaternions  (vectors) 
is  not  a  pure  quaternion  (vector)  since  the  scalar  part  of  the  product  is  usually 
not  zero.  Sometimes  the  product  of  two  quaternions  with  non-zero  scalar  part 
yields  a  pure  quaternion.  For  example,  from  (4.26) 


e  e  = 


> 

=  ~  l 


e  e 


(4.28) 


This  is  only  possible  if  the  above  product  is  hyperimaginary,  i.e.,  a  pure 
quaternion  or  vector. 

We  now  try  to  formulate  a  rotation  operator  in  terms  of  quaternions  operating 
on  a  vector  (i.e.,  a  quaternion  with  a  zero  scalar  part) .  The  simplest  thing  to  try 
is  multiplication  of  a  vector  from  the  right  or  left  by  a  unit  quaternion.  A  unit 
quaternion  is  chosen  since  rotations  preserve  the  length  or  norm  of  a  vector,  and 
we  can  see  from  (4.21)  that  multiplication  by  a  unit  quaternion  will  not  change 
length.  First  we  will  try  multiplication  from  the  left.  We  shall  see  that 
multiplication  from  the  left  only  (or  from  the  right  only)  is  unsatisfactory. 


We  choose  for  the  rotation  operator  the  unit  quaternion 


f  -  I  f  -1 

X  =  ^  XQ  ,  X  j  =  [  cos(p)  ,  sin(P)  X  j 


(4.29) 


where  the  “hat”  denotes  a  unit  vector  and 

2  -  -  2  2  2  2 
XQ  +  X<*"qx.  ‘0  +  +  X2  +  S 

For  the  vector  we  choose  a  quaternion  representation  with  zero  scalar  part 


(4.30) 


q  = 


o  .  ?] 


(4.31) 


Then  we  try  representing  a  rotation  by 


For  quaternion  q'  to  be  a  “vector”,  the  scalar  part  must  vanish.  But  multiplying 
a  pure  vector  by  a  quaternion  produces  another  quaternion  with  a_non-vector 
component.  Thus,  unless  we  make  an  orthogonality  assumption  that  k.q  —  0  so 
the  scalar  part  vanishes,  the  quaternion  multiplication  does  too  much.  But  such 
an  assumption  is  too  restrictive.  Alternatively,  we  could  try  multiplying  by  X 
from  the  right.  This  gives 

q’  =  qX  =  [o  ,  4]  [x0  ,  x]  =  [~?.X  ,  XQ4  +  q  X  X  ]  (4.33) 

Note  that,  from  (4.22)  we  can  write 

q'=  qX'1  =  [+?.X,-X0q  +  XX9"  ]  (4.34) 

The  scalar  parts  of  (4.32)  and  (4.34)  have  opposite  sign.  This  suggests  we  might 
try  combining  (4.32)  and  (4.33)  into  a  similarity  transformation  in  the  hope  that 
we  might  be  able  to  get  the  scalar  part  to  drop  out.  This  strategy  turns  out  to  be 
a  good  one.  We  shall  see  that  this  approach  does  not  require  any  restrictive 
assumptions  such  as  the  orthogonality  of  the  quaternion  and  the  vector.  In 
addition,  it  will  turn  out  to  be  equivalent  to  the  matrix  rotation  operator 


described  in  (3.5). 

q'=  Xq  X*1 

(4.35) 

Expanding 

11 

> 

O 

>*\ 

» _ i 

O 

i— j 

o 

-x  3  = 

(4.36) 

r-  2_ 

\^\.q  +  kQq»\  +  XX^.X  ,  +  XQ<?  +  XflXX^  +  \.qk  —  \QqX\  —  XX^xXj 

Using  the  identities  [A  x  B  J.  A  —  0  and  /\xBxC  =  B[A.C^  —  with  the 

normalization  condition  for  the  unit  quaternion  X,  viz.,  (4.29),  this  becomes 


Thus  q'  defined  by  (4.35)  is  still  a  vector  or  pure  quaternion  and  its  length  is 
preserved,  as  required  for  a  rotation.  We  need  to  manipulate  these  terms  so  that 
q  appears  on  the  right  with  some  operator  expression  to  its  left.  This  vector  part 
can  be  written  in  matrix  component  form  by  using  the  following  identities 


r-  _1 

LXX«J. 


e. .,  X  .q. 

ijk  j  k 


(4.38) 


Also, 


o 

1 

>* 

OJ 

X2 

X3  0 

-X2  X, 

0 

1^ 

III 

L  JX.<?  JX  J  -  \k 

i 

Jr 

II 

X1X1  X1X2 

X1X3 

X2Xl  X2X2 

X2X3 

X3Xl  X,X2 

X3X3 

^  2 
«3 


\  ( - T  )  _ 


«2 

«3 


(4.39) 


where  the  superscript  denotes  the  transpose  (interchange  of  rows  and 
columns).  Thus  we  can  write 


-  (  2  )  t  - "  _ 

q'_  =  i(_2Xa-lJ/  +  2  A.  A..  +  2X0Aj5 


(4.40) 


=  2 


A0  +  Ax  -  % 


XiX2  +  X0X. 


^  j  ^  ^  Xi  Q  ^ 


^1^2  ^0^3 

x0  +  x2  -  ^ 


x2x3  +  x0x1 


^2^5  Vi 


x0  +  x3  -  % 
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=  2 


^Xp  +  Xj'-xJ-Xj2] 
X  ^  X  qX  ^ 


XjX3  ^0^2 


XjX.  x0x. 


%[x0‘-x;+x2-x;] 

^2^3  "*"  ^0^1 


X|X ^  4*  X.q^** 


X2X, 


%[x0-x;-x:+x3] 


-  T  q 


where  I  is  the  identity  or  unit  matrix.  The  second  version  of  the  matrix  was 
obtained  using  (4.29).  The  expression  (4.35)  which  involves  quaternion 
operations  is  equivalent  to  (4.40)  which  involves  matrix  operations. 

Alternatively,  the  quaternion  nature  of  (4.35)  may  be  used  more  directly  as  a 
4x4  matrix  by  using  the  two  expressions  for  quaternion  multiplication  given  in 
(4.10)  and  making  use  of  the  normalization  condition  (4.19). 


q'  =  X  q  X'1  =  X  q  X  = 


(4.41) 


+  xo 

"X2 

-X3 

+  Xo 

+  X1 

+  X2 

+  X 

+X1 

+  xo 

"X3 

+  X2 

"X1 

+  Xo 

"X3 

4*  X 

+  X2 

+  X3 

+  xo 

"X2 

+  x3 

+  xo 

-X 

+  X3 

"X2 

+  X1 

+  Xo 

"X3 

~X2 

+  XI 

+  x 

1 

3 

l2 

ll 

l0 


0 


\q 


*2 

*3 


1 

0 

0 

0 

0 

0 

2  2  2  2 
[X0  +  X1  _X2-X3] 

2[X1X2  "  X0X3] 

2  2  2  2 

2fXlX3  +  W 

0 

2[XxX2  +  X0X3) 

[X0-Xl+X2-X3^ 

2[X2X3  -  X^j) 

2  2-2  2 

^  2 

0 

2[X2X3  -  XQX2] 

2[X2X3  +  X0Xj] 

(Xo"x1  —  ^2  +  x3) 

«3 

How  do  we  see  that  (4.40)  or  (4.41)  correspond  to  (3.5)?  Let  us  start  with 
(4.29). 
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(4.29) 


X  =  cosO)  ,  sin(P)  n 


where 

XQ  =  cos(0) 

1/2 

f  2  1 

sin(p)  =  l  1  ~  x0  J 
ik1  +  ;X2  +  kk3 

n  = 

sin(p) 

We  will  establish  the  connection  between  (4.35),  (4.40)  and  (3.5)  by  means  of  an 
example  rather  than  a  comprehensive  and  rigorous  proof,  by  dealing  with  a 
rotation  about  a  single  one  of  the  coordinate  axes.  Of  course,  a  complete  general 
rotation  can  be  built  up  by  several  such  component  rotations.  Take  the  special 
case  where  n  lies  along  the  y-axis.  Then  XQ  =  cos(p),  X2  =  sin(fi)  and  Xj  —  X3 
=  0. 

The  reader  can  readily  verify  by  substitution  into  (4.40)  that  (3.3)  will  be 
generated  except  for  the  curious  fact  that  the  angle  <J>  is  equal  to  2|3!  This 
provides  insight  into  the  interpretation  of  the  rotation  associated  with  a 
quaternion.  An  Euler  quaternion  generates  a  rotation  about  an  axis  determined 
by  the  vector  part.  The  half  angle  of  rotation  is  determined  by  the  arc  tangent  of 
the  ratio  of  the  length  of  the  vector  part  to  the  scalar  part. 

The  matrix  T  is  the  rotation  matrix  from  moving  axes  to  inertial  axes.  It  is 
numerically  identical  to  the  much  less  symmetric  and  computationally  more 
complex  expression  in  (3.5).  But  as  simple  as  this  quaternion  form  of  the 
rotation  matrix  appears  to  be,  recall  it  is  just  the  matrix  and  vector  parts  of 
(4.37)  under  matrix  algebra.  But  (4.37)  is  identical  to  (4.35),  which  is  the 
expression  of  a  rotation  of  a  vector  represented  as  a  quaternion  with  zero  scalar 
part.  It  is  nothing  more  than  a  unit  quaternion  X  multiplied  by  the  vector 
multiplied  by  the  inverse  of  the  unit  quaternion  X,  where  multiplication  means 
quaternion  multiplication.  Examine  (3.5)  and  note  how  cumbersome,  and 
computationally  awkward  it  is.  Then  examine  (4.40)  and  finally  (4.35).  How 
deceptively  simple,  elegant  and  computationally  efficient  (4.40)  and  (4.35)  are! 
By  now  the  reader  should  begin  to  have  an  appreciation  for  the  power  of  the 
quaternion  formalism.  There  is  a  draw  back,  however.  One  can  easily 
intuitively  grasp  the  Euler  angles.  The  four  quaternion  components  are  not  so 
easily  subject  to  intuition. 

Because  of  this,  it  is  more  convenient  to  define  the  initial  conditions  of  a 
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simulation  in  terms  of  Euler  angles  rather  than  the  quaternions  themselves.  We 
can  use  the  Euler  angles  to  generate  a  rotation  matrix  T.  It  then  remains  to  find 
the  quaternions  that  would  generate  the  same  rotation  matrix.  We  already  know 
how  to  generate  T  from  the  quaternions  X..  How  do  we  do  the  inversion,  obtain 
the  X.  from  T? 

I 

The  expression  (4.35)  can  be  inverted  to  give  the  quaternions  in  terms  of  the 
rotation  matrix  T  in  (4.40).  This  inversion  can  be  accomplished  as  follows. 
From  the  diagonal  elements  of  (4.40)  and  the  unit  quaternion  normalization 
condition  (4.29),  we  obtain  the  following  four  simultaneous  equations,  where  we 
will  be  using  the  trace,  Tr  [  T  ],  which  is  defined  as  the  sum  of  the  diagonal 
elements  of  the  matrix  T,  viz.,  Tr  [  T  ]  =  T  n+T  2z+T  33' 


+ 


\  + 
2 

2 

X1  + 


(4.42) 


These  can  be  solved  simultaneously  to  give  the  following. 


1/2 

x0  =  ±%[  1  +  7>(T)) 

1/2 

Xj  =  ±%[  1  +  2TU  -  7>(T)] 

1/2 

X2  =  ±fc[  1  +  IT 22  -  7>(T)] 

1/2 

X3  =  ±M  1  +  27 33  -  Tr(T)] 


There  is  a  sign  ambiguity  to  resolve.  The  chief  constraint  on  the  algebraic  signs 
of  the  X  is  that  the  rotation  matrix  T  from  which  the  X  were  generated  should  be 
regenerated  when  these  X  are  substituted  in  (4.40).  Note  that  no  negative  matrix 
elements  can  arise  if  all  the  X  have  the  same  sign,  whether  positive  or  negative. 
Thus  the  choice  of  sign  is  not  trivial.  The  diagonal  elements  of  T  pose  no 
constraint  since  the  X  occur  only  as  squares  in  the  diagonal.  The  off-diagonal 
elements  appear  as  cross  terms  and  are  the  key  to  our  task.  From  (4.40),  we  can 
take  all  possible  combinations  of  T„±  T„,  i^j. 
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\v- 


*1  Tn  + 


r2.1 


VS  “  ~  r12> 

=  %t  t13  +  r31] 

*0^2  =  y*\T  j3  -  ^3^ 

X2X3  ~  ¥*[T 23  + 


(4.44) 


Vl  %^~32  ^23^ 

In  generating  these  expressions  to  investgate  the  sign  ambiguity,  we  see  that 
(4.44)  provides  an  alternate  form  for  determining  the  X,  providing  that  one  of 
the  kj  is  known.  For  example,  suppose  we  obtain  Xfl  from  (4.43).  Then  XJ?  X 
and  X3  can  be  obtained  from  (4.44)  by  dividing  the  appropriate  expressions  in 
(4.44)  by  XQ.  Thus  we  obtain 


XQ  =  ±%  1  +  Tr[T]  ^ 


4X, 


T  -  T 

32  23 


X„  = 


—  f 


4X, 


[T13 


T 
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X.  = 


4X, 


21 


-  T 


12  J 


(4.45) 


This  hybrid  solution  combining  (4.43)  and  (4.44)  has  an  unexpected  advantage. 
The  sign  ambiguity  remains  only  for  XQ.  This  sign  ambiguity  is  not  significant 
since  changing  the  sign  of  XQ  will  change  the  sign  of  all  the  X  ..  But  the 
quaternions  always  appear  as  products  of  pairs  of  quaternions  or  the  Square  of  a 
quaternion.  Thus  the  rotation  matrix’is  invariant  tinder  a  sign  change  for  X  , 
See  (4.40).  Although  the  .sign  convention  doesn’t  matter,  it  is  conventional 
chosen  to  be  positive. 

The  above  inversion  can  have  numerical  problems  in  computation,  however. 
The  X  .  for  i=  1,2,3,  become  ill-defined  when  the  kQ  in  the  denominator  vanishes. 
This  Will  occur  if  the  trace  of  the  rotation  matrix  equals  -1,  which  would  happen 
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when  i|j  =  0  and  4>  =  -ir  radians  or  180  degrees.  Furthermore,  numerical  round 
off  errors  occur  in  the  immediate  neighborhood  of  this  region  when  using  (4.45) 
for  computation. 

It  is  not  difficult  to  avoid  this  problem  since  we  did  not  have  to  start  with  AQ  in 
(4.43).  Another  choice  could  have  been  made  before  invoking  (4.44)  for  the 
other  3  A..  This  could  be  done  by  evaluating  all  four  A.  in  (4.43),  picking  the 
largest  and  then  using  this  dominant  A.  to  find  the  other  A  in  (4.44).  We  would 
again  obtain  results  that  would  be  insensitive  to  sign  ambiguity  but  would  now  be 
insensitive  to  computational  round  off  errors  as  well  by  putting  us  far  away  from 
any  singularity. 

This  algorithm  is  described  in  Table  1  with  two  modifications.  For  consistency, 
the  signs  are  examined  at  the  end  of  the  algorithm.  If  A0  is  negative,  the  signs  of 
alF  A.  are  reversed  to  keep  AQ  positive,  according  to  our  convention.  We  have 
alrea'dy  seen  that  changing  the  overall  sign  of  all  four  quaternions  has  no  effect 
on  the  rotation  matrix.  Secondly,  evaluating  all  four  quaternions  using  (4.43)  is 
computationally  wasteful.  It  is  sufficient  to  compare  the  trace  and  the  three 
diagonal  elements  of  T  to  select  the  dominant  quaternion.  The  complete 
algorithm  is  shown  in  Table  1. 


In  order  to  use  the  quaternions  we  must  have  a  scheme  for  determining  the  time 
rate  of  change  of  the  quaternions  so  that  these  rates  can  be  integrated  to  provide 
updated  quaternions  as  the  system  evolves  in  time.  These  would  be  expressions 
that  play  a  role  completely  analogous  to  (3.16)  through  (3.18)  for  the  Euler 
angles.  This  derivation  could  have  been  accomplished  by  standard  matrix 
algebraic  techniques.  However,  the  derivation  would  have  been  exceedingly  long 
and  tedious.  Using  quaternions,  it  is  rather  simple  and  straight  forward. 

Since  q  and  q'  are  vectors  or  pure  quaternions,  they  obey  the  transformation 
law  given  in  (4.35), 

-l  .  -l  ,  . 

q'  =  A  q  A  q=AqA 

(4.46) 

q  =  A  q  A  q  —  A  q  A 


Now  we  wish  to  take  the  time  derivative  of  q  in  the  moving  coordinate  frame.  As 
is  well-known,  a  term  must  be  added  due  to  the  angular  velocity  II  of  the 
rotating  coordinate  frame.  Taking  the  derivative 
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Table  1.  Obtaining  the  Quaternions  from  the  Rotation  Matrix  T 


•  Define  T00  =  Tr  [  T  1. 

•  Compare  T„  where  1  =  0,1, 2,3. 

Find  the  dominant  quaternion  (most  positive  or  least  negative.) 

The  index  of  the  dominant  T{.  determines  the  dominant  X.  to  be  used  in  (4.43). 

•  Determine  the  other  three  X^ ,  j  #  i  in  (4.44). 

The  four  cases,  with  the  dominant  X  first  are 


4A... 


*2  = 
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1  4*  Tr  (T  )] 
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k3  =  +%[  1  +  2T3J  +  Tr(T)]’'2 
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4x. 


f  T  V  7 
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Xl'-' 
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K+t„] 


l  r  '  v 

1  r  i 

X3  =  7“  T23+T« 

X2=~T23+r32 

4X,  1 

4A.  1  1 

•  Examine  the  algebraic  sign  of  XQ.  If  negative,  change  the  sign  of  all  four 
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d q  d  _i  -  _ 

-  =  — [  X  q'  X]  +  [  0  ,  ftx$] 

dt  dt 

=  X  2q'  X  +  X  1q  'X  +  X  V  X  +  [  0  ,  ftx<?]  (4.47) 

All  these  terms  are  pure  quaternions  or  vectors.  Since 

X_1  q  '  X  =  X_1  X  q  X~*  X  =  q  (4.48) 

we  conclude  that 

X  *q'  X  +  X_1q'X  =  — (0,  ftxq  ]  =  +  [  0  ,  qxfi  ] 

or  (4.49) 

X  XXqX_1X  +  X  *XqX  *X  =  +[  0  ,  qxft  ] 

.  -  1  —  1  * 

=  X  Xq  +  qX  X 

- 1  -  ’  ■  - 1  - 1  - 

Differentiating  X  X  =  [1,0]  gives  X  X  =  —  X  X.  Thus  from  (4.7) 

[  0  ,  qxft  ]  =  q  X  1  X  —  X  q  =  2  [  0  ,  qx[X_1  X]  ]  (4.50) 

-l  •  . 

Recall  from  (4.28)  that  X  X  is  a  pure  quaternion  or  vector.  Thus 

1  -l  • 

-  ft  =  X  X  (4.51) 

2 

We  use  the  i,  j,  k  notation  described  in  (4.1)  through  (4.6),  since  this  formalism 
makes  it  easier  to  group  terms  and  find  a  matrix  operator  equivalent.  We  obtain 

2X_1X  =  ft  - 

[  0  ,  ft  ]  =  [  0  +  iftj  +  jn2  +■  kft3  ]  =  (4.52) 

2  [X0-/X1-;X2-*X3]  [XQ  +  *X1  +  ji2  +  k\3) 
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2  [VoniV  ^2^2+  ^3^3  J 


+  2i  |—  XaX.0+ XqXjH- X3X2~ ^2^3]  + 
+  2/:[-\3X0  +  X2\1-X1X2+X0X3] 


This  may  be  expressed  into  the  matrix  form  (recall  equation  4.10) 


0 

+  Xo 

+  X1 

+  x2 

+  X3 

xo 

ft 

X 

=  2 

"Xi 

+  xo 

+  x3 

X2 

X1 

ft 

y 

-X2 

"X3 

+  xo 

+  X1 

X2 

ft 

Z 

"X3 

+  X2 

”X1 

+  X0 

X3 

(4.53) 


The  top  element  in  the  column  matrix  for  the  angular  velocity  can  be  seen  to 
vanish  by  taking  the  time  derivative  of  the  normalization  condition  on  X,  viz., 
(4.29).  Similarly,  the  first  line  of  (4.52)  can  be  quaternion-multiplied  from  the 
right  by  X  to  yield 


X  =  %Xft  = 


%[X0  +  iX1  +  ;X2  +  kX3][0  +  iftx  +  j&y  +  kClJ  = 


(4.54) 


%  [  yj-x^-x^-x^]  +  WIXjO+X^-X^  +  X^  ] 

+  %j  [  X20+X3ftjc  +  Xony-X1Oz]  +  [X30— Xjfl^  +  Xjfl^  +  ] 


This  may  be  organized  into  the  matrix  form  using  (4.10) 


+  xo 

"X1 

"X2 

-X„ 

0 

0 

+  X1 

+  xo 

"X3 

+  x2 

... 

ft 

X 

+  X2 

+  x3 

+  Xo 

X 1 

ft 

y 

+  X3 

"X2 

+  X1 

+  xo 

ft 

1 

(4.55) 
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See  equation  (4.10)  Note  that  both  (4.53)  and  (4.55)  are  free  of  singularities, 
unlike  their  analogs  (3.19)  and  (3.20).  Both  matrices  (4.53)  and  (4.55)  are  equal 
to  a  unit  matrix  times  X  added  to  an  antisymmetric  matrix.  These  matrices  are 
inverses  of  one  another  as  well  as  transposes  of  one  another.  Hence  they  are 
orthogonal  and  .preserve  the  length  of  the  vectors  they  operate  on.  Thus  the 
vectors  (  \1  X  X2  X3  )  and  (  0  XI  £1^  )  are  related  by  a  rotation  in  4-space 
with  the  latter  having  only  a  three-vector  zpart.  Compare  these  elegant  properties 
and  the  simplicity  and  low  computational  burden  associated  with  the  use  of  these 
equations  versus  (3.19)  and  (3.20).  If  the  above  derivation  is  not  simple  enough 
for  you,  the  matrix  in  (4.55)  may  be  obtained  by  realizing  it  is  the  inverse  of  the 
matrix  in  (4.53),  and  can  be  obtained  by  substituting  into  the  matrix  in  (4.53)  the 
inverse  of  the  unit  quaternion  X.  But  from  (4.22)  we  realize  that  the  inverse  of  a 
unit  quaternion  is  obtained  merely  by  changing  the  sign  of  the  last  three  ( 
“vector-like”)  components  of  the  quaternion.  Thus  all  that  needs  to  be  done  to 
invert  this  matrix  is  to  change  the  sign  of  the  off-diagonal  elements.  There  is 
something  here  for  everyone:  the  mathematician,  the  physicist,  the  engineer  and 
the  programmer. 

Since  quaternions  are  not  as  intuitive  as  Euler  angles,  it  is  sometimes  desirable  to 
move  back  and  forth  between  the  Euler  angle  and  quaternion  representations. 
Going  from  Euler  angle  to  quaternions  representation  can  be  achieved  by  using 
the  Euler  angles  to  evaluate  the  transformation  matrix  using  equation  (3.5).  The 
transformation  matrix  T  can  then  be  put  into  (4.41)  to  obtain  the  quaternions. 
Note  that  after  (4.41),  it  was  noted  that  the  sign  for  Xo  is  not  important.  Why 
this  is  so  becomes  apparent  shortly.  From  the  first  row  and  last  column  of  (3.5), 
and  (4.40) 


—  tt/2  <  0  <  tt/2 


—  7!  <  «|>  <  TT  (4.56) 


0  <  4>  <  2tt 


Note  that  theta  can  only  be  defined  in  this  way  over  a- range  it  where  the  arcsin  is 
defined,  viz.  between  4- ir/2  and  —tt/2.  However,  and  <f>  can  be  found  over  a 
full  360  degrees  or  2  rr  radians  by  taking  into  account  the  sign  of  the  numerator 
and  the  sign  of  the  denominator  of  the  last  two  expressions.  See  the  algorithm 
programmed  in  Table  2.  Note  further  that  either  sign  selection  for  \o  made  in 
(4.41)  wii]  yield  back  the  same  Euler  angles  in  the  above  expression.  Thus  the 
sign  choice  is  arbitrary. 


sin0  =  -731  =  2[X0X2- XjXj] 


T21  2^1X2+X0X3^ 


tarn]/  = 


li 


2  2  2  2 

N,  +  ^1  ~~  ^2  3 
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tancj)  = 


33 


2[x2X3+x0x1] 


2  2  2  2 
X0~X1_X2  +  X3 
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Table  2.  Evaluation  of  ARCTAN  (A  ,B)  Over  All  Four  Quadrants 


IF  B>0 
IF  B=  0,  A>0 
IF  B=  0,  A<0 
IF  B<0 


tan  *(A/B) 


it/2 


—  tt/2 

it  +  tan  (A/B) 


The  denominator  in  two  of  the  expressions  given  in  (4.56)  can  vanish  when  or 
6  are  at  ±ir/2.  The  algorithm  in  Table  2  will  treat  this  correctly. 

We  shall  see  in  the  next  chapter  that  plane-fixed  coordinates  will  require  4>  =  0. 
When  applied  to  the  expression  for  tab  4>  in  (4.56),  this  constraint  will  further 

require 


X2k3  + 


Vl  =  0 


(4.57) 


A  note  of  caution  should  be  addressed,  viz.,  the  necessity  for  renormalization. 
Round  off  error  in  a  computer  simulation  gradually  causes  the  normalization 
condition  (4.20)  or  (4.29)  to  fail.  Thus  the  quaternions  need  to  be  regularly 
renormalized  in  a  simulation. 
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5  EQUATIONS  OF  MOTION 


In  this  section,  the  rigid  body  equations  of  motion  will  be  developed  for  three 
types  of  coordinate  systems:  body-fixed,  aeroballistic  and  plane-fixed.  Body- 
fixed  coordinates  rotate  (roll,  pitch,  yaw)  with  the  body  of  a  projectile.  Hence 
the  angular  velocity  ft  of  the  coordinate  axes  is  equal  to  the  angular  velocity  to  of 
the  projectile  body.  See  (3.13)  through  (3.15).  For  body-fixed,  coordinates 

ft  =0)  =  P  =  —  ijr  sin  0  +  <j>  (5.1) 

XX 

ft  =  co  =  Q  =  sin  <i>  cos  0-1-0  cos  <j> 
y  y 

ft  =o)  =  R  =  \ b  cos  <b  cos  0  —  0  sin  4* 

Z  1 

Recall  that  in  aerodynamics  the  components  of  the  angular  velocity  of  the 
projectile  body  o>  are  conventionally  denoted  by  P,  Q  and  R.  These  coordinates 
are  the  natural  choice  for  guided  projectiles  since  the  seeker  and  sensor  outputs, 
actuator  parameters,  and,  so  forth  are  most  simply  and  naturally  expressed  in 
body  coordinates.  However,  for  spin  stabilized  projectiles,  a  very  small 
integration  time  step  is  required.  Otherwise  the  coordinate  system  will  roll 
through  too  great  an  angle  during  the  time  step  and  smear  the  direction  of  forces 
such  as  gravity.  To  deal  with  this,  plane-fixed  coordinates  are  utilized.  These 
coordinates  pitch  and  yaw  with  the  projectile  but  do  not  roll  with  it.  In 
particular,  one  axis  is  constrained  to  remaining  in  a  single  plane.  In  our  case,  the 
y-axis  will  be  constrained  to  the  horizontal  plane.  See  equations  (3.22)  and 
following.  This  implies  that  the  roll  Euler  angle  of  the  plane -fixed  frame  satisfies 
the  relations 


d$/dt  =  0  ,  4>  =  0  (5.2) 

Thus,  although  the  projectile  is  rolling,  the  fixed  plane  coordinate  system  is  not 
rolling,  in  the  sense  that  <}>  vanishes  for  the  fixed  plane  axes.  However,  the  x 
component  of  the  angular  velocity  of  the  frame,  viz.  ft^  does  not  vanish  and  ft^ 
of  the  frame  does  not  equal  P  (the  x  component  of  the  projectile  angular 
velocity).  These  relations  may  be  found  in  equations  (3.26)  through  (3.28). 

ft  =  —  vh  sin0  =  —R  tan0 

X 

Cly <2-  =  'i  (5.3) 

ft  =  R  =  cos0 

Z 

Comparing  this  expression  for  ft  with  (5.1),  we  see  that  4>  is  the  roll  rate  of  the 
projectile  with  respect  to  pfane-fixed  coordinate  frame.  (5.1)  carT  be 
reconstructed  from  (5.3)  by  adding  <|>  to  (5.3)  and  rotating  with  (3.3).  This 
expression  in  useful  if  the  Euler  angle  representation  is  used.  With  the 
quaternion  representation  for  plane-fixed  coordinates,  we  make  use  of  (3.5), 
(4.40)  and  4>  =  0  to  obtain 
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—  2  [x.X,  -  XX  1 
—  T’  I  1  3  0  2| 

31 


tan0  = 


33 


2  2  2  2 
X0  “  X1  “  X2  +  x3 


(5.4) 


Thus 


ft  = 


2 R  [x^  -  XQX2 


2  2  2  2 
X0  -  X1  -  X2  +  S 


(5.5) 


The  denominators  in  (5.4)  and  (5.5)  can  vanish,  corresponding  to  0  =  ±ir/2.  In 
defining  plane-fixed  coordinates,  we  had  to  impose  Euler  angle  type  algebra  on 
the  quaternions  and  have  corrupted  them  with  Euler  angle  type  singularities. 
Thus,  quaternions  in  plane-fixed,  coordinates  should  not  be  used  in  vertical 
trajectories.  Use  body-fixed  coordinates  instead.  This  is  not  burdensome 
computationally  since  the  direction  of  gravity  is  along  the  axis  of  roll.  Finally, 
we  could  choose  ft^  —  0  in  (^5.1)  for  aeroballistic  coordinates.  This  choice  will 
simplify  equations  of  motion 

In  the  following  development,  results  will  be  derived  using  the  term  ft^ .  This 
will  allow  all  three  coordinate  formalisms  to  be  developed  simultaneously.  At  the 
end,  the  results  can  be  specialized  to  one  or  the  other  by  letting  ft^  equal  P  for 
the  body-fixed  case  or  equal  zero  for  the  aeroballistic  ("zero  P")  case  or  equal 
—  R  tan0  (or  its  quaternion  equivalent  -  see  (5.3)  and  (5.5)).  ft  =  Q  and  ft  = 
R  for  all  three  cases.  From  (3.13)  we  see  that  the  definitions*  for  aeroballistic 
and  body-fixed  frames  are  not  equivalent  but  become  the  same  in  the  limit  of 
small  0  or  small  dty/dt .  The  advantage  of  choosing  plane-fixed  coordinates  for  a 
non-rolling  system  is  that  it  has  no  y  component  of  gravity  in  a  flat  earth  model, 
thus  eliminating  the  possibility  of  gravity  smearing  due  to  roll  during  integration 
of  the  equations  of  motion.  (This  will  be  shown  in  equation  (5.10)  below.) 


In  summary,  the  body-fixed  coordinates  roll,  pitch  and  yaw  with  the  projectile 
and  act  as  if  physically  attached  to  the  projectile.  Plane-fixed  coordinates  pitch 
and  yaw  with  the  projectile  but  do  not  roll  with  it.  The  y-axis  is  constrained  to 
move  in  the  horizontal  plane.  See  the  discussion  for  equations  (3.22)  to  (3.33). 
The  Euler  angle  rotation  matrix  for  the  plane-fixed  case  can  be  obtained  from 
the  body-fixed  matrix  (3.5)  by  letting  <j>  =  0.  Equivalently,  (3.3)  is  replaced  by 
the  unit  identity  matrix. 


*  Vaughn,  Harold  R.,  "A  detailed  Development  of  the  Tricyclic  Theory,"  Sandia  Laboratories,  SC-M-67-2933, 
Albuquerque,  NM,  1968. 
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The  form  of  Newton’s  law  F  =  dmv/dt  is  valid  only  in  inertial  (i.e.,  non¬ 
accelerating)  coordinate  frames.  If  the  coordinate  system  is  rotating,  Newton’s 
law  will  not  be  valid  because  a  rotation  is  an  acceleration.  However,  the  law  can 
be  amended  for  rotational  frames.  As  is  well-known,  Newton’s  law  for  linear 
accelerations  and  forces  in  a  rotating  frame  takes  the  form  (the  superscript  M 
denotes  the  moving  frame,  body-fixed,  plane-fixed  or  aeroballistic)1234 

F  +  mg  =  mV  +  £lxmV  —  mV  +  milxV  (5.6) 

dt 

F  contains  applied  forces  such  as  thrust  and  aerodynamic  forces.  Since  the 
coordinate  system  is  non-inertial,  it  also  contains  "fictitious”  terms  such  as 
centrifugal  and  Coriolis  "forces".  Derivatives  of  inertial  properties  such  as  mass 
will  be  omitted  in  this  development1 * * 4 5.  Defining  V *  =  U,  V  =  V,  and  =  W, 
the  components  of  the  above  vector  equation  is  now 

F^  +  mg*  =  mU  +  m[QW  —  RV] 

F  +mgM  =  mV  +  m[RU -D.  W]  (5.7) 

y  y  x 

F '  4-  mgM  —  mW  4-  m  V  —  QU  ] 

Rearranging 

t  * 

U  =  - +  g*-QW+RV 

m 


1  John  H.  Blakelock,  "Automatic  Control  of  Aircraft  and  Missiles,**  John  Wiley  and  Sons,  New  York,  1965. 

^  Keith  K.  Symon,  "Mechanics, *'  Addison -Wesley,  Reading,  Mass,  1960. 

Goldstein,  Herbert.  "Classical  Mechanics",  p  136,  Addison  Wesley.  Reading,  Mass,  1959. 

4  Landau,  L.  D.;  and  E.  M.  Lifshitz,  "Classical  Mechanics",  p  128,  Addison  Wesley,  Mass,  1960. 

5  The  thrust  of  a  reaction  engine  that  is  measured  in  a  test  stand  already  contains  the  effects  of  the  rate  of  change  of  the 
mass.  Since  this  information  is  usually  available  for  input  into  a  simulation  rather  than  nozzle  pressures,  derivatives  of 
the  mass  do  not  appear  in  the  equations  of  motion.  However,  if  the  thrust  is  to  be -reconstructed  from  pressure 
measurements  at  the  nozzle,  mass  derivative  terms  and  the  velocity  of  the  exhaust  gases  would  have  to  be  taken  into 
account  in  the  equations  of  motion.  This  will  be  discussed  in  detail  in  a  future  report. 
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F 

v  =  —+gM-Ru  +  a  w 

m 

Fz 

W  =  —  +  g^-SlxV  +  QU 

m 


(5.8) 


The  components  for  g  can  be  obtained  by  multiplying  the  gravity  vector  in  the 
earth  frame  by  the  appropriate  T  matrix.  For  the  special  case  for  a  flat  earth,  g 
only  has  a  vertical  or  z  component  pointing  downward.  (Coriolis  and  centripetal 
acceleration  corrections  should  be  made  since  a  flat  earth  is  not  really  inertial 
since  the  earth  rotates.  If  distances  flown  and  time  of  flight  are  short,  these 
corrections  are  negligible.) 

We  need  the  inertial  to  body  transformation  (  from  the  inverse  of  (4.40)  )  and 
the  representation  of  the  gravity  vector  in  a  flat  earth  as  (  0  ,  0  ,  g  ) 


•2[A.2X3+X0X1)g  for  flat  earth  (5.9) 

2  2  2  2 
tXo“Xi  _X2  +  X31'? 

If  the  user  wishes,  the  components  for  g  in  the  aeroballistic,  plane-fixed  or 

body-fixed  frames  for  a  flat  earth  could  be  obtained  in  terms  of  the  Euler  angles 

-IT. 

instead  of  quaternions  by  using  the  expression  for  T  =  T  in  (3.5)  instead  of 
(4.40). 

M  .  „ 

gx  =  -g  sine 

gM  =  +£sin<}>  cos0  for  flat  earth  (5.10) 

M 

g  =  +  g  cosq)  COS0 


M 
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Note  that  there  is  no  plane-fixed  y-axis  component  of  gravity  since  sin  <j>  for 
plane-fixed  coordinates,  advantage  of  the  definition  of  plane-fixed  coordinates 
adopted  in  this  development.  It  would  not  be  true  if  =  0  was  chosen.  In 
general,  if  g  has  x  or  y  components  in  inertial  coordinates,  the  full  rotation 
matrix  T  would  have  to  be  used.  See  (3.5)  and  (4.40). 

Recalling  that  O  is  the  angular  velocity  of  the  coordinate  frame  with  respect  to 
the  inertial  frame,  to  is  the  angular  velocity  of  the  body  itself  with  respect  to  the 
inertial  frame,  and  the  moment  of  inertia  tensor  is  symmetric,  i.e.,  /  = 

Newton’s  laws  for  angular  velocities  and  moments  are  of  the  form 


M  =  L  +  SlxL  =  /  to  +  ftjt[/co)  (5.11) 
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Time  derivative  of  the  angular  velocity  of  the  body 
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cross  product  Angular  velocity  of  the  body 
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In  aerodynamic  conventional  notation. 


M  =  L 


M  =  M 
y 

M  =  N 

i 


(5.12) 


Body-fixed,  aeroballistic  and  plane-fixed  coordinates  can  be  developed 
simultaneously  by  writing  ft^  =  P  for  body-fixed,  ft^  =  0  for  aeroballistic  and 
by  writing  for  plane-fixed  coordinates 


a  =  2R[\1x3-x0x2)/[x02-\12-x22+x32]  (5.13) 


—  —  R  tan6 


See  (5.2)  or  (5.5). 
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(5.14) 
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Q  [-1  P  -I  Q  +/  R  ]  -  R  [-1  P  +1  Q  -I  R  ) 

zx  zy  zz  yx  yy  yz 

R  [  +1  P  -1  Q  -1  R  ]  -  n  [  -/  P  -I  Q  +1  R  ] 

xx  xy  xz  x  zx  zy  zz 

O  [  -/  P  +r  Q  -I  R  ]  -  Q  [  +/  P  -1  Q  -I  R  ] 

x  yx  yy  yz  xx  xy  xz 


This  may  be  written 


L  =  I  P 

XX 


(5.15) 


[/  -I  1 
l  yy  J 


QR  <  -  Vanishes  with  axial  symmetry 


+  Iyz  [RZ-Q2\  +  Izz  [-<2/>-«j+/t>  <-  Vanish  i: 


Vanish  if  products  of  inertia  vanish 


M  =  /  PR  +  I  Q  -  I  R£l 

xx  yy  zz  x 


+  /  [ga  -/?]+/  fa  p-r2\+i  -qr-p 

yz  I  x  J  xz  I  x  J  xy  ^  J 


<-  Vanish  if  products  of  inertia  vanish 


N  =  -I  PQ  +  /  n  Q  +  I  R 

xx  yy  x  zz 


+,„  [-n„*-e]+/„  W-A*',,  W~^p] 


<-  Vanish  if  products  of  inertia  vanish 
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If  fi  is  replaced  by  P,  equations  (5.14)  or  (5.15)  are  the  general  equations  of 
motion  for  the  angular  velocity  components  of  a  rigid  body  in  body-fixed 
coordinates.  For  convenience  in  programming,  terms  that  vanish  when  the 
products  of  inertia  vanish  have  been  grouped  together.  Body-fixed  coordinates 
are  appropriate  for  unsymmetric  bodies  and  for  guided  projectiles,  since  the 
sensors  and  control  system  are  naturally  described  in  body-fixed  coordinates  that 
roll  with  the  body.  Note  that  the  components  of  the  moment  of  inertia  tensor  I 
will  not  change  due  to  the  rotation  of  the  body-fixed  coordinates. 


For  axially-symm etric  spin-stabilized  projectiles,  plane-fixed  coordinates  are 
more  appropriate.  If  body-fixed  coordinates  are  used  with  rapidly  spinning  spin 
stabilized  projectiles,  the  integration  time  step  in  a  6  degree-of-freedom 
simulation  is  driven  to  be  very  small,  increasing  the  simulation  run  time.  This  is 
necessary  to  keep  the  roll  angle  during  the  integration  .time  step  small.  This 
avoids  smearing  gravity  in  the  expressions  for  x,  y  and  z.  See  (5.7)  or  (5.8). 
With  plane-fixed  coordinates  and  axial  symmetry,  the  products  of  inertia  7  ,  7 , 
and  7  vanish,  and  7  =  7  .  is  replaced  using  (5.13)  in  (5.14)  or  (5.15). 

Note  tfiat  the  components  of“he  moment  of  inertia  tensor  I  would  generally  vary 
with  time  if  the  projectile  were  rotating  but  the  frame  were  not.  Not  only  would 
this  be  a  complication  but  it  would  drive  down  the  integration  time  step  and 
increase  integration  time.  Similarly  for  the  aerodynamics.  Since  the  motivation 
for  plane-fixed  coordinates  is  greater  computational  speed,  it  is  pointless  to 
eliminate  gravitational  smear  and  substitute  inertial  or  aerodynamic  smear.  Thus 
axial  symmetry  in  mass  properties  and  aerodynamics  is  generally  assumed. 


After  one  more  result  is  obtained,  these  formulae  will  be  collected  in  Tables  3  to 
8.  From  (4.55),  using  the  appropriate  expression  for  we  can  write  down  the 
time  derivatives  of  the  quaternions  for  these  frames. 
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For  body-fixed, 

\  =  as  -  *s] 

X,  =  %[+FX0  -  Q\3  +  /?X2] 
X2  =  %[+PX3  +  Q\q  -  /?xj 

x3  =  %[-px2  +  exx  +  /?x0] 

For  plane-fixed,  using  (4.55)  and  (5.5) 


(5.16) 


X0  = 


2R\, 


X1X3  X0X2  1 


r  2  2  2  2,  . 
[X0~Xl~X2  +  X3l 


-  Q\^  ■- 


•) 


x3  * 


x,e  + 


[  2  2  2  2_ 
X0-X1-XI  +  X!]j  j 


\x  =  14 


+  2.RXr 


XjXj  XQX. 


^xJ-xJ-Xj+Xj2] 


-  ex 


3  +  ^2j 


(5.17) 


=  -% 


X.  J? 


x.o  + 


[x0!-x;-x2:+x;]j  _ 
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X,  = 


+  2  R\, 


X1X3  X0^2 


2  2, 


[Xq  Xj  Xj  +  X.] 


+  Q\q  -  R\1 


=  +% 


Xj  R 


Xn£?  + 


2  2  2  2, 1 
x0-Xi-x2  +  x3]j 


1 

r 

^1^3  ^0^2 

'i 

>" 

II 

-2*X, 

2  t  2  2 

+  exj  +  *x0 

1 

l 

*  [X0  ~~  ^1  ~  ^2  +  ^3 1  " 

j 

=  +% 


xo* 


x,G  + 


[xJ-xf-x-j+xJl]  J 


A  singularity  will  occur  if  the  denominator  vanishes.  From  (3.22)  and  (4.40), 
this  denominator  is  just  T  =  cos(0),  which  vanishes  for  0  =  ±tt/2.  This  is  the 
same  singularity  discussed  after  (5.5).  It  is  not  advisable  to  use  plane-fixed 
coordinates  for  near  vertical  trajectories.  Use  body-fixed  coordinates  instead. 
For  vertical  trajectories,  gravity  smearing  due  to  roll  rate  is  not  the  problem  it  is 
for  other  trajectories. 

All  the  quaternion  and  Euler  angle  results  for  the  linear  and  angular  equations  of 
motion  are  collected  in  the  following  tables  for  body-fixed  and  for  plane-fixed 
coordinates  respectively. 

For  completeness,  we  note  that  the  aeroballistic  frame  equations  of  motion  can 
be  obtained  by  letting  ft  vanish  in  (5.8)  and  (5.15),  and  letting  ft  vanish  in 
(4.55)  or  P  vanish  in  (*5.16).  As  with  plane-fixed  coordinates,  we  further 
simplify  by  neglecting  the  products  of  inertia.  This  form  will  be  found  in  Table 
5.  Further  simplification  results  from  assuming  axial  symmetry  and  constant 
mass.  See  Tables  5  and  6. 
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Table  3.  Body-Fixed  Equations 


M 

U  =  - +  g  -QW  +RV 

m 

F 

y  M 

V  =  - +  g  -RU+PW 

°y 

m 

F 


(5.1),  (5.8) 


W 


M 


+  g  -PV  +  QU 


m 


-l 


The  components  of  gravitational  acceleration  are  obtained  using  T 
For  the  flat  earth  approximation,  use  (5.9)  and  (5.10). 


L  =  I  P 

XX 


[,  -/  1 

L  «  yy  J 


QR  <  -  Vanishes  for  axial  symmetry 


(5.15) 


+  /  R2-Q2\  +  1  -QP-R  +/  RP-Q 

y2  L  J  xi  «.  J  xy  L  J 


<-  Vanish  if  products  of  inertia  vanish 


M  —  +1  Q+'I  -/  1 RP 


yy 


L «  « J 


r  i  r  2  2i  f  1 

+  I  \QP  “  R  +  7  [P  -R  J+/  —  QR  —  P  J  <-  Vanish  if  products  of  inerti 


inertia  vanish 


W  =  */  -I  PQ  +  /  R 

yy  xx  zz 


r  1  r  1  2  21 

+  /  —  PR  —  Q  +/  QR  —  P  j  +  /  Q  —P  I  <-  Vanish  if  products  of  inertia  vanish 

yz  J  xzi  JxyL  J 
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Table  4.  Time  Development  of  the  Body-Fixed  Transformation 

Matrix  Parameters 


X0  =  %["PX1  "  “  *X3j 

=  %  [+/>X0  "  Qk3  +  R  X2  ] 
X2  =  %[+PX3  +  Q\  -  /?X2] 

X3  =  %[-PX2  +  QX1  +  /?X0] 


i 

4> 

8 


or 

[Q  sin( <f>)  +  R  cos (<j>)  J 
cos (9) 

P  +  [Q  sin(fy)  +  R  cos(<|>))  tan{ 8) 
Q  c os  ( 4> )  -  R  sin(<$>) 


(5.16) 


(3.16) 

(3.17) 

(3.18) 


( The  expressions  (3.16)  and  (3.17)  have  a  singularity  at  ±  tt/2.  This  singularity 
does  not  occur  in  the  quaternion  expression  (5.16).  See  discussion  in  text.) 
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Table  5.  Plane-Fixed  Equations 
(Axial  Symmetry,  Iyy  =  Izz=  /  ,  products  of  inertia  vanish) 

F 

x  M 

U  =  - +  g  -QW+RV 

X 

m 

F 

_ y_  M 

V  =  - +  g  -RU+CIW  (5.8) 

m 

F 

.  ^ 

W  =  - +  £  -ft  V  +  0t/ 

z  x 

m 

The  components  of  gravitational  acceleration  are  obtained  using  T  \ 

For  the  flat  earth  approximation,  use  (5.9)  and  (5.10). 

L  =  I  P 

XX 

M  =  ItQ  +  I^PR  -  1R&X  (5.15) 

N  =  1  R  -  I  QP  +  /  ft  Q 

t  X X  14  t  X  ^ 


where  f l  is  obtained  from 

X 

ft  —  —  R  tan8 

jr 

or 


ft 

X 


2  R 


[x,*3  - 


\>V 


9 


(5.3) 


(5.5) 
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Table  6.  Time  Development  of  Plane-Fixed  Transformation 

Matrix  Parameters 


-%  M2  + 


[  2  2  2  2  1 
xo-xi~x:  +  Mj 


M2  + 


(5.17) 


[\02-Xi2-\:2+X32]]  J 


+  %  + 


[xJ-x’-x’+X,2]]  j 


+%  I  Xjg  + 


[x0:-x;-x;+x;i]  j 


<}>  =  0 

<t>  =  o 


cos (6) 


e  =  q 


(3.16) 


(3.18) 


(The  above  expressions  (5.17)  and  (3.16)  are  singular  near  the  vertical.  Use  body- 
fixed  coordinates  instead.  See  the  discussion  in  the  text.) 
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Table  7.  Aeroballistic  Equations 


These  are  obtained  by  letting  XI  =  0,  and  assuming  no  products  of  inertia.  A 
further  simplification  can  be  made  by  letting  /  =  =  /f  (axial  symmetry). 


F 

.  X  w 

U  =  +  gx  -QW+RV 

m 

F 

y 

v  =  +  g“-RU  (5.1),  (5.8) 

m 

F 

z  M 

W  =  +  gz  +QU 

m 


The  components  of  gravitational  acceleration  are  obtained  using  T 
For  the  flat  earth  approximation,  use  (5.9)  and  (5.10). 


L  =  lj>  (5.15) 

M  =  +1  Q+I  RP 

t  XX 

N  =  +//?-/  PQ 

t  XX 
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Table  8.  Time  Development  of  the  Aeroballistic  Transformation 

Matrix  Parameters 


K  =  »[- 

Xj  =  nf  -  qx3  +  «x2] 
x2“  *[  +  ex0  -  RXj] 
x3  =  *[  +  ex,  +  «x0j 


or 

4» 

4> 

0 


(i2  sin  (<j>)  +  R  coj(4>)) 
cos (0) 

{Q  5zn(4>)  4-  R  cos (<j>))  tan(Q) 
»J/sin(0) 

Q  cos(<$>)  —  R  sin(<$>) 


(The  Euler  angle  expressions  (3.16)  and  (3.17)  have  a  singularity  at 
Equation  (5.16B)  comes  from  letting  P=0  in  (5.16). 


(5.16B) 

(3.16) 

(3. 17), (5. 1) 

(3.18) 

ir/2.) 
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6.  INTEGRATION  OF  EQUATIONS  OF  MOTION 


6.1  Plane-Fixed  Equations 

Recall  for  the  plane-fixed  case,  the  force  equations  from  Table  5  are 


F 

X  M 

U  =  - +  gx  -QW  +RV 

m 

F 

y  M 

V  =  +  g  —  RU  +  ft  W 

y  x 

m 

F 

_ f_  M 

W  =  - +  gx  -ilxV  +  QU 

m 


(6.1) 


where  is  given  by  (5.13)  and  gM  js  given  by  (5.9)  for  a  flat  earth  or  more 

M  -1  / 

generally  by  g  =  T  g  ,  where  the  subscript  I  refers  to  an  inertial  frame. 
These  expressions  are  readily  integrated  numerically. 


For  the  plane-fixed  case,  the  moment  equations  from  Table  5  can  be  put  into  an 
uncoupled  form  for  integration,  where/  =  /  =  /  . 

t  yy  zz 

L 

P  =  - 

I 

XX 


Q  =  Im-I'RP+I  Rn-  1 

I  XX  t  X  I 

I  .  J 


R  = 


N+I  QP 


IQ  & 

t  X 


These  are  readily  integrated  numerically. 


(6.2) 
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6.2  Body-Fixed  Equations 


For  body-fixed  case,  ft  =  P  and  from  (5.8) 


U  =  —  +  *"-  QW+RV 
m 


y  M 

V  =  - +  g  -RU  +  PW 

m 

F 

_ f_  B 

W  =  - +  gz  -PV  +  QU 

m 


which  are  readily  integrated  numerically. 


For  the  body-fixed  case,  equations  (5.15)  may  be  written 


-  ~ [i+V2 +/„*-/, {r.Q,x}] 


1 


where 


\ , 


■r»  o  »n 


/ 

xv 

s  -  —  [w+/„p+/,;e-/3(/>.e.R)J 

I 

1Z 


f  =  \l  -/  }qR+1  [i?2—  <22 ]  —  /  QP+I  RP 

J  1  [  zz  yy  J  yz  <-  J  zx  xy 


(6.3) 


(6.4) 
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(6.5) 


/,  =  [/  -7  ]/?P+7  QP+I  \p2-R2\-J  QR 

2  [  xx  zz  J  yz  ^  zx  ^  ->  xy  ^ 

/  =  [/  -/  jpg-/  PP+7  0P+7  f(22-P2| 

3  [  yy  xx  J  yz  zx  ^  xy  J 


As  written,  these  equations  are  coupled.  If  the  products  of  inertia  vanish,  the 
numerical  integration  is  straight  forward  because  the  equations  become 
uncoupled  in  the  derivatives,  as  in  (6.2)  above.  Even  if  they  do  not  vanish,  the 
products  of  inertia  7^  ,  I  ^  and  7  are  generally  quite  small.  This  suggests  a 
simple  approximation*  Tfie  equations  could  be  solved  by  using  the  derivatives 
dP/dt,  dQ/dt  and  dR/dt  on  the  right  side  of  equation  (6.2)  from  the  previous  time 
step.  Since  the  products  of  inertia  are  typically  small,  this  approximation  should 
be  adequate  in  practice.  For  more  precision,  the  results  could  be  iterated.  That 
is,  the  results  for  the  derivatives  dP/dt,  dQ/dt  and  dR/dt  on  the  left  side  could  be 
put  back  into  right  side  to  obtain  a  better  approximation  before  performing  an 
integration. 


If  the  products  of  inertia  are  not  small  or  if  one  wishes  to  avoid  this 
approximation,  equations  (6.4)  and  (6.5)  can  be  put  into  the  form 


L-f.  =  +7  P-7  Q-I  R 

1  xx  xy  ^  xz 

M-f  =  -I  P+1  Q-I  R  (6.6) 

2  x y  yy  yz  v  ' 

N  — /  =  -7  P-7  Q+I  R 

3  xz  yz  zz 


They  can  be  solved  simultaneously  to  uncouple  the  derivatives  of  P,Q  and  R  by 
inverting  the  matrix  of  moment  of  inertia  components.  Formally  we  can  write 
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-1 


p 

+  1 

XX 

-/ 

xy 

-i 

xz 

L~f ,  1 
1  1 

Q 

= 

-1 

xy 

+ 1 

yy 

-1 

yz 

s: 

i 

to 

R 

-1 

XZ 

-/ 

yz 

+  / 

zz 

N-f,  | 

Evaluating  this  inverse  is  somewhat  tedious  but  the  procedures  for  inverting  a 
matrix  are  well  known  *.  If  we  denote  the  matrix  to  be  inverted  by  J 


J"1  =  [  det  J  r1 


I  I  -I 

yy  zz  yz 

i  i  +  i  / 

xz  yz  xy  zz 

i  i  +i  i 

xy  yz  xz  yy 


11+11 

xz  yz  xy  zz 

i  i  -i2 

XX  zz  xz 

11+11 

xy  xz  xx  yz 


11+11 

xy  yz  xz  yy 

11+11 

xy  xz  xx  yz 

i  i  -i2 

xx  yy  xy 


(6.8) 


where 


det  J  =  111  -Ill  -III  -  I  l2  -  I  I2  -  I  I2 

xx  yy  zz  xy  yz  xz  xz  xy  yz  yy  xz  zz  xy  xx  yz 


(6.9) 


*  Gelb,  Arthur,  et  al.,  Applied  Optimal  Estimation  Theory",  p  17,  MTT  Press,  Cambridge,  Mass,  1974. 
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6.3  Aeroballistic  Equations 


Recall  for  the  aeroballistic  case,  the  force  equations  from  Table  7  are 


F 

x_  M 

U  =  +  gx  —  QW  +  RV 

m 

F 

y  M 

V  =  - +  g  -RU 

y 

m 

F 

_ z  M 

W  =  —  +  g T  +QU 
m 


(6.10) 


The  moment  equations  from  Table  7  are 


L 

P  —  — 
1 

XX 


Q  -  \m -I  Rp] 

j  L  "  J 

!/V  +/  QP 

l  XX 


(6.11) 


These  are  readily  integrated  numerically. 
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APPENDIX 


ALGORITHMS  FOR  IMPLEMENTATION 
OF  THE  EQUATIONS  OF  MOTION 
IN  SIX  DEGREE  OF  FREEDOM  COMPUTER  SIMULATIONS 
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TRANSITIONS  BETWEEN  EULER  ANGLES  AND  QUATERNIONS 

The  Initialization  Problem 


Since  most  people  are  more  intuitively  comfortable  with  Euler  angles  than  with 
quaternions,  the  Euler  angle  to  quaternion  transformation  can  be  used  to  input 
initial  conditions  in  Euler  angle  format  for  the  convenience  of  the  user  and 
convert  to  quaternions  for  internal  use  in  a  simulation  if  so  desired.  Conversely, 
quaternions  used  internally  by  a  simulation  can  be  converted  to  Euler  angles 
prior  to  generating  output,  for  the  convenience  of  the  user. 


QUATERNIONS  TO  EULER  ANGLES: 

The  Euler  angles  can  be  evaluated  directly  from  the  quaternions  or  indirectly 
from  a  rotation  matrix  that  had  been  developed  from  the  quaternions.  Use 
(4.56). 

• 

Note  that  the  denominator  of  the  expressions  for  4*  and  for  0  in  (4.56)  6an 
vanish.  The  algorithm  in  Table  2  in  this  document  can  handle  this  case  correctly. 


EULER  ANGLES  TO  QUATERNIONS: 

1)  Evaluate  the  transformation  matrix  T  from  the  Euler  angles  using  (3.5). 
(With  plane-fixed  coordinates,  the  roll  Euler  angle  <f>  must  be  set  to  zero: 
Alternatively,  (3.22)  can  be  used.) 

2)  Evaluate  the  quaternions  using  Table  1. 
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TRANSITIONS  BETWEEN  BODY-FIXED,  PLANE-FIXED 
AND  AEROBALLISTIC  COORDINATES 


Plane-fixed  coordinates  are  more  efficient  for  modeling  spin-stabilized, 
unguided,  rotationally  symmetric  projectiles.  There  is  no  component  of  gravity 
outside  of  the  x-z  plane  in  a  flat  earth  model.  Thus  these  coordinates  are 
insensitive  to  gravity  smearing  because  of  roll.  However,  plane-fixed 
coordinates  have  a  singularity  for  vertical  trajectories  and  body-fixed  coordinates 
are  preferred  for  such  trajectories.  Body-fixed  coordinates  are  also  more 
appropriate  for  guided  stages  or  other  stages  that  don’t  have  the  required 
symmetry.  Similarly,  aeroballistic  coordinates  are  also  more  efficient  for  axially- 
symmetric,  spin-stabilized  projectiles  than  body-fixed  coordinates.  Furthermore, 
the  equations  of  motion  are  simpler  for  this  choice  than  the  other  two  candidate 
coordinate  frames. 

This  document  permits  the  development  of  a  6  degree  of  freedom  simulation  in 
which  the  coordinate  frame  can  be  changed  from  one  stage  to  another  to  use  the 
coordinate  frame  that  is  most  appropriate  or  efficient  in  each  particular  stage  of  a 
trajectory  simulation.  Generally,  other  than  changing  the  equations  of  motion, 
nothing  special  needs  to  be  done  when  transitioning  between  one  type  of  frame 
and  another.  There  is  an  exception  for  transitioning  to  plane-fixed  coordinates 
from  other  frames. 

When  transitioning  to  plane-fixed  coordinates  there  is  a  discontinuous 
change  in  the  orientation  of  the  coordinate  system,  since  4>  must  vanish  in 
plane-fixed  coordinates.  This  requires  the  following  adjustments: 

1)  The  projectile  angular  velocity  vector  (P,Q,R)  must  be  temporarily 
transformed  to  non-moving  (i.e.,  inertial)  coordinates  using  the  last  value 
of  the  rotation  matrix. 

2)  A  new  rotation  matrix  T  must  be  generated.  (If  the  quaternion 
representation  is  being  used,  the  equivalent  Euler  angles  must  be 
regenerated  first,  as  shown  in  (4.56).)  Set  coordinate  frame  Euler  roll 
angle  <j>  to  zero,  retaining  the  regenerated  pitch  and  yaw  angles. 
Recalculate  the  rotation  matrix  with  these  new  Euler  angles  using  (3.5)  or 
(3.22). 

3)  Use  this  matrix  to  rotate  the  projectile  angular  velocity  vector  in  the 
non-moving  frame  back  to  the  moving  (plane-fixed)  frame  to  obtain  the 
new  P,  Q  and  R. 

4)  If  using  the  quaternion  formalism,  regenerate  the  quaternions  from  the 
rotation  matrix  using  Table  1. 

Resume  calculations  using  the  appropriate  equations  of  motion  for  the  type  of 
coordinate  frame  being  used,  for  either  the  Euler  angle  or  quaternion 
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representations  described  elsewhere  in  this  document.  The  pitch  0  and  yaw  ijj 
Euler  angles  for  the  body  and  for  plane-fixed  axes  will  be  identical.  The  plane- 
fixed  frame  roll  Euler  angle  <f>  is  a  constant  and  identically  zero.  The  roll  Euler 
angle  for  the  aeroballistic  frame  is  not  constant  but  generally  will  not  vary  much 
from  its  value  when  transition  occurred. 

When  using  plane-fixed  or  aeroballistic  coordinates,  the  body  angular  velocity 
component  P  must  be  reconstructed  for  use  with  the  aerodynamics.  When 
transitioning  back  to  body-fixed  coordinates,  a  roll  angle  can  be  restored  by  first 
constructing  the  inverse  (i.e.,  transpose)  of  the  roll  rotation  matrix  from  (3.3). 
The  new  full  roation  matrix  is  obtained  by  taking  the  matrix  product  of  the 
existing  plane-fixed  rotation  matrix  and  the  roll  matrix,  in  that  order.  This 
procedure  will  not  corrupt  the  other  variables  in  the  equations  of  motion. 
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BODY-FIXED  COORDINATES  USING  QUATERNIONS 


Rotation  matrix  from  body-fixed  to  inertial  coordinates: 

Calculate  T(  kQ,  X^  X2>  X3)  from  (4.40). 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body-fixed,  is  obtained 
by  taking  the  transpose  of  this  matrix. 


Time  derivatives  of  the  quaternions: 

Calculate  \Q,  Xr  X2,  S  from  (5.16). 


Other  Equations  of  Motion  (See  Tables  3  and  4.) 

a)  Force  equations:  Use  (5.8)  with  fl*  =  P. 

See  (5.1).  Obtain  components  of  gravity,  in  body-fixed  frame  from  the  inertial 
frame  by  using  T  .For  flat  earth,  use  (5.9). 

b)  Moment  equations:  Use  (6.4)  and  (6.5)  or  (6.9)  to  (6.11). 
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PLANE-FIXED  COORDINATES  USING  QUATERNIONS 


Rotation  matrix  from  plane-fixed  to  inertial  coordinates: 


Calculate  T(  XQ,  Xr  X2,  X3  )  using  (4.40) 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body-fixed,  is  obtained 
by  taking  the  transpose  of  this  matrix. 


Time  derivatives  of  the  quaternions: 

Calculate  X0,  X  ,  X9,  X3  from  (5.17). 
with 

ft  « - 

x  2  2  2  2 

\  ~  X1  ~  X2  +  X3 

from  (5.5).  There  is  a  constraint  in  (4.57). 


Other  Equations  of  Motion  (See  Tables  5  and  6): 


a)  Force  equations:  Use  (5.8)  .with  ft^  from  (5.5).  See  above.  - 

See  (5.1).  Obtain  components  of  gravity  in  body-fixed  frame  by  using  T  \  For 
flat  earth,  use  (5.9).  Because  of  the  constraint  (4.57),  there  is  no  y  component 
of  gravity. 

b)  Moment  equations:  Use  (6.2). 
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AEROBALLISTIC  COORDINATES  USING  QUATERNIONS 


Rotation  matrix  from  body-fixed  to  inertial  coordinates: 

Calculate  T(  kQ,  ky  kr  kj  from  (4.40). 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body-fixed, 
by  taking  the  transpose  of  this  matrix. 


Time  derivatives  of  the  quaternions: 

Calculate  iQ,  X,,  kr  X3  from  (5.16B),  or  (5.16)  with  P=0. 


Other  Equations  of  Motion  (See  Tables  7  and  8.) 

a)  Force  equations:  Use  (5.8)  with  11^  =  P. 

See  (5.1).  Obtain  components  of  gravity  in  body-fixed  frame  from 
frame  by  using  T  \  For  flat  earth,  use  (5.9). 

b)  Moment  equations:  Use  (6.11). 


is  obtained 


the  inertial 
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BODY-FIXED  COORDINATES  USING  EULER  ANGLES 


Rotation  matrix  from  plane-fixed  to  inertial  coordinates: 

Calculate  T(  <J>,  0,  ^  )  from  (3.5). 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body-fixed,  is  obtained 
by  talcing  the  transpose  of  this  matrix. 


Time  derivatives  of  the  Euler  angles: 


ft  =  p 

X 

ft  =  Q  ft  =  R 

y  z 

(5.1) 

(q  sin(<f>)  +  R  cos(4>)j 

(3.16) 

* 

co.r(0) 

<i> 

P  +  [q  sin(4>)  4-  R  cos (<$>))  tan(0) 

(3.17) 

0 

Q  co5(<j>)  —  R  sin(4>) 

(3.18) 

Other  Equations  of  Motion  (See  Tables  3  and  4.) 

a)  Force  equations:  Use  (5.8)  with  ft^  =  P. 

See  (5.1).  Obtain  components  of  gravity  in  body-fixed  frame  by  using  T  .  For 
flat  earth,  use  (5.10). 

b)  Moment  equations:  Use  (6.4)  and  (6.5)  or  (6.7)-(6.9). 
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PLANE-FIXED  COORDINATES  USING  EULER  ANGLES 


Rotation  matrix  from  plane-fixed  to  inertial  coordinates: 

Calculate  T(  <J>=0,  0,  i}»)  from  (3.5). 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body- fixed,  is  obtained 
by  taking  the  transpose  of  this  matrix.  Note  that  <j>  is  zero  in  (3.5). 


Time  derivatives  of  the  Euler  angles: 

Use  <{>  =  <j>  =  0  instead  (3.17)  and  =  Q  and  =  R. 

R 

ij*  =  from  (3.16). 

cos (8) 

0  =  Q  from  (3.18) 


Other  Equations  of  Motion  (See  Tables  5  and  6.) 

a)  Force  equations:  Use  (5.8)  with  from  (5.3). 

See  (5.1).  Obtain  components  of  gravity  in  plane-fixed  frame  by  using  T  \  For 
flat  earth,  use  (5.10)  with  4>  =0. 

b)  Moment  equations:  Use  (6.2). 
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AEROBALLISTIC  COORDINATES  USING  EULER  ANGLES 


Rotation  matrix  from  plane-fixed  to  inertial  coordinates: 

Calculate  T(  cf>,  0,  )  from  (3.5). 

The  inverse  rotation  matrix,  from  inertial  coordinates  to  body-fixed,  is  obtained 
by  taking  the  transpose  of  this  matrix. 

(5.1) 

(3.16) 

(3.17) 

(3.18) 

Other  Equations  of  Motion  (See  Tables  7  and  8.) 

a)  Force  equations:  Use  (5.8)  with  =  0. 

See  (5.1).  Obtain  components  of  gravity  in  body-fixed -frame  by  using  T  \  For 
flat  earth,  use  (5.10). 

b)  Moment  equations:  Use  (6.11). 
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CONSTRAINTS,  LIMITATIONS  AND  PRACTICAL  REQUIREMENTS 

I.  Normalization  constraint  on  quaternions: 

We  require  XQ  +  X1  +  X2  +  X3  —  1  from  (4.29). 

At  each  integration  time  step  (or  at  least  at  frequent  intervals), 
normalize  by  dividing  each  X.  by 

-  1/2 

N  =  [xfl2+  X2+  X2+  X2] 


II.  Constraint  for  plane-fixed  coordinates: 


We  require  X2  X3  +  ~  ®  from  (4.57). 

Check  this  constraint  regularly.  If  it  begins  to  fail: 

a)  regenerate  the  Euler  angles  from  the  matrix  T  using  (4.56), 

b)  set  <j>  =  0,  and 

c.)  regenerate  quaternions  from  (4.42). 


III.  Euler  angle  singularity: 


Terminate  simulation  if  6  is  too  close  to  ±90  degrees  because  of  the  singularity 
at  that  angle.  See  (3.16)  and  (3.17).  (N.B.,  the  Euler  angle  rotations  of  roll, 
pitch  and  yaw  may  be  chosen  to  move  the  singularity  to  occur  along  the 
horizontal  rather  than  the  vertical  axes.  A  better  solution  is  to  use  quaternions 
with  body-fixed  coordinates.  For  plane-fixed  coordinates,  see  paragraph  IV 
below. 


IV.  Quaternion  singularity  for  plane-fixed  coordinates: 


This  singularity  is  similar  to  III  except  it  only  occurs  for  plane-fixed  coordinates 
and  not  for  body-fixed  or  aeroballistic  coordinates.  Body-fixed  coordinates 
chould  be  used  for  vertical  trajectories  rather  than  plane-fixed. 
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V.  Axial  symmetry  requirement  for  plane-fixed  and  aeroballistic 
coordinates: 


a)  /  =  / 

yy  « 

b)  Products  of  inertia  7^,  /„,  7yz>  lyx,  7„.  7Z),  aI1  vamsh‘ 

c)  Aerodynamic  coefficients  do  not  depend  on  roll  angle. 
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